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Abstract 

The  precision  navigation  capabilities  of  the  Global  Positioning  System  (GPS)  are 
used  extensively  within  US  military  operations.  However,  GPS  is  highly  vulnerable 
to  intentional  and  unintentional  external  interference.  Therefore,  a  need  exists  to 
develop  a  non-GPS  precision  navigation  method  to  operate  in  GPS  degraded  envi¬ 
ronments.  This  research  effort  presents  the  theoretical  limits  of  a  precision  navigation 
method  based  on  an  inertial  navigation  system  (INS)  aided  by  angle  measurements 
with  respect  to  lunar  surface  features  observed  by  a  fixed  camera.  To  accomplish  this 
task,  an  extended  Kalman  filter  (EKF)  was  implemented  to  estimate  INS  drift  errors 
and  bring  in  simulated  lunar  feature  angle  measurements  to  correct  error  estimates. 
The  research  scope  focused  solely  on  the  feasibility  of  lunar  vision  aided  navigation 
with  INS  where  only  measurement  noise  effects  from  a  simulated  CCD  camera  and 
barometer  were  considered.  Various  scenarios  based  on  camera  specifications,  lunar 
feature  quantity,  INS  grade,  and  lunar  orbital  parameters  were  conducted  to  observe 
the  INS  drift  correction  by  lunar  feature  angle  measurements.  The  resulting  trade 
spaces  presented  by  the  scenarios  showed  theoretical  substantial  improvement  in  the 
navigation  solution  with  respect  to  a  stand  alone  INS. 
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THEORETICAL  LIMITS  OF 


LUNAR  VISION  AIDED  NAVIGATION 
WITH  INERTIAL  NAVIGATION  SYSTEM 

I.  Introduction 


1.1  Purpose 

The  development  of  non-Global  Positioning  System  (GPS)  based  precision  navi¬ 
gation  methods  is  vitally  important.  This  is  due  to  GPS’s  high  vulnerability  to  inten¬ 
tional  and  unintentional  external  interference  coupled  with  the  US  military’s  heavy 
reliance  on  GPS.  The  following  research  develops  the  theoretical  limits  of  a  navigation 
method  integrating  inertial  sensors  and  lunar  feature  angle  measurements.  The  goal 
is  to  provide  immutable,  robust,  self-contained,  autonomous  precision  navigation  in 
GPS  degraded  environments. 

Since  ancient  times,  the  ability  to  travel  from  one’s  present  location  to  an  intended 
destination  has  been  practiced  on  land  and  sea  [28].  Over  time,  this  ability  and  its 
associated  methods  became  the  field  of  study  called  navigation  derived  from  the 
Latin  word  navis  for  sea  faring  ship  [3]  [12].  Of  the  ancient  navigation  methods 
employed,  around  two  millennia  ago,  the  Polynesians  are  known  to  have  developed 
a  form  of  celestial  navigation  whereby  determining  one’s  position  based  on  careful 
observation  of  celestial  bodies  [28].  The  Chinese  and  the  Italians  are  recorded  to 
have  independently  developed  the  compass  during  the  12th  and  13th  centuries  CE 
[28]  [34],  Later  developments  in  celestial  navigation  techniques  in  Europe  during  the 
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18th  century  allowed  for  latitude  and  longitude  to  be  estimated  with  the  development 
of  the  sextant  and  the  chronometer  respectively  [34] . 

Technological  advances  around  the  beginning  of  the  20th  century  brought  about 
new  possibilities  for  precision  navigation  methods.  New  understanding  behind  ro¬ 
tational  motion  and  acceleration  paved  the  way  for  navigation  based  upon  inertial 
measurements  in  the  form  of  an  inertial  navigation  system  (INS)  [28].  The  manipu¬ 
lation  of  radio  waves  in  the  1930s  proved  that  navigation  by  radio  waves  was  another 
viable  method  of  navigation  known  as  the  American  long  range  navigation  system 
(Loran)  [9]. 

With  the  advent  of  the  space  age  in  the  1950s,  the  US  and  Soviet  Union  realized 
the  capability  to  navigate  by  their  own  manufactured  celestial  bodies  (or  satellites) 
rather  than  rely  on  natural  celestial  bodies  (e.g.  moon  and  stars).  The  navigation 
satellite  constellation  developed  by  the  US  was  known  as  Navstar  GPS.  The  first  GPS 
satellite  was  launched  in  1978  [9] .  The  core  positioning  methodology  of  GPS  is  based 
upon  streaming  and  receiving  transmitted  GPS  signals  from  the  satellite  constellation 
to  ground  receivers  [12].  The  detected  timing  differences  due  to  the  signal  propagation 
from  satellite  transmitter  to  a  ground  based  receiver  can  be  multiplied  by  the  speed 
of  light  to  get  a  pseudorange  [12].  A  pseudorange  is  a  range  containing  error  from  the 
ground  based  receiver’s  position  to  the  satellite  transmitter’s  position  [12].  A  position 
fix  can  be  attained  from  knowing  the  receiver’s  respective  timing  of  received  signals 
in  relation  to  at  least  four  satellites  [9]. 

The  Gulf  War  in  1991  acted  as  the  debut  of  GPS  in  an  armed  conflict  through 
the  use  of  GPS  based  precision  guided  munitions  and  ground  unit  navigation  [25]. 
GPS’s  level  of  precision  and  robust  navigation  qualities  are  said  by  McLendon  to 
have  enabled  the  “left  hook”  maneuver  by  Coalition  forces  which  destroyed  the  Iraqi 
armored  divisions  [22] .  By  the  end  of  the  20th  century,  GPS  had  become  an  integral 
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navigation  system  for  the  US  military  [27].  GPS  was  not  without  fault  and  weakness. 
Over  time,  the  system  has  been  found  to  be  degraded  in  certain  circumstances.  Urban 
canyons  and  jamming  systems  are  known  to  subdue  the  effect  of  GPS  in  said  areas 
[9]  [12].  Thus,  the  need  to  counter  such  weakness  exists  and  can  be  achieved  through 
the  development  of  other  navigation  means  independent  of  GPS. 

1.2  Problem  Statement 

This  research  effort  assumes  GPS  has  been  degraded  and  other  means  of  trajectory 
estimation  are  required  for  an  aircraft  to  properly  navigate.  Trajectory  is  defined  as 
the  position,  velocity,  and  attitude  of  a  vehicle  over  time  [30]. 

The  non-GPS  precision  navigation  system  chosen  for  this  research  effort  is  an  INS. 
An  INS  performs  dead  reckoning  navigation  meaning  that  a  trajectory  estimate  can 
be  made  based  on  an  INS’s  initial  position  in  combination  with  direction  and  speed 
measurements  [28].  Direction  and  speed  measurements  are  given  by  accelerometers 
and  gyroscopes  which  are  inertial  sensors  within  an  INS.  Accelerometers  measure  spe¬ 
cific  forces  acting  upon  a  vehicle  body  and  gyroscopes  measure  changes  in  rotational 
velocity  also  with  respect  to  a  vehicle  body  [28].  INSs  are  highly  accurate  for  short 
time  periods  but  accumulate  drift  error  for  an  extended  time  length  [5].  An  optical 
sensor  with  a  view  of  the  Moon  can  provide  aiding  angular  measurements  to  decrease 
the  accumulated  drift  error  in  the  INS. 

The  goal  of  this  research  effort  is  to  explore  the  theoretical  improvement  upon 
an  INS  estimated  trajectory  aided  by  angular  measurements  with  respect  to  lunar 
features  observed  by  a  camera.  The  improvement  of  the  trajectory  estimation  will 
be  measured  in  terms  of  the  reduction  of  position  state  errors  over  time.  The  air¬ 
craft  trajectory,  INS  errors,  lunar  features,  and  optical  sensor  measurements  will  be 
simulated. 
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1.3  Research  Contributions 


The  following  are  the  research  contributions  for  this  thesis: 

1.  Model  the  drift  errors  for  an  INS  (commercial,  tactical,  and  navigation  grade). 

2.  Construct  a  model  of  an  Earth-based  aircraft  trajectory  with  features  simulated 
at  lunar  surface  locations. 

3.  Simulate  INS  drift  errors  aided  by  lunar  feature  angle  measurements  to  demon¬ 
strate  a  decrease  in  trajectory  error  over  various  scenarios  where  measurement 
noise  effects  are  only  considered  (e.g.  INS  grade,  pixel  size,  lens  focal  length, 
feature  quantity,  and  lunar  distance  and  geometry). 

As  only  measurement  noise  effects  upon  the  system  are  considered  within  this  research 
effort,  the  results  in  Chapter  4  are  not  indicative  of  realistic  lunar  vision  aided  navi¬ 
gation  with  INS.  However,  the  results  present  plausible  trade  spaces  for  lunar  vision 
aided  navigation  with  INS  that  exist  with  respect  to  the  various  scenario  parameters 
mentioned  above  in  the  research  contributions. 

1.4  Research  Assumptions 

To  limit  the  scope  of  the  research  effort  to  consider  only  measurement  noise  effects 
upon  lunar  vision  aided  navigation  with  INS,  the  following  simplifying  assumptions 
were  made: 

1.  The  lunar  feature  angular  measurements  observed  by  the  camera  were  mod¬ 
eled  with  unbiased  Gaussian  white  noise.  This  research  effort  does  not  include 
the  more  complicated  error  contained  in  more  realistic  feature  extraction  and 
matching  techniques  such  as  the  Scale  Invariant  Feature  Transform  algorithm 
(SIFT). 
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2.  In  a  realistic  scenario,  there  will  be  error  in  the  camera  alignment  with  respect 
to  the  local  navigation  frame.  However,  in  this  research  effort,  it  is  assumed 
that  the  alignment  to  the  horizon  and/or  fixed  stars  was  available  and  known. 
Therefore,  the  camera  alignment  error  was  not  included. 

3.  The  camera  was  modeled  as  fixed  to  the  vehicle  body.  For  a  dynamic  vehicle, 
this  scenario  would  typically  induce  feature  smearing  over  the  duration  of  the 
camera  exposure.  However,  this  affect  was  not  considered  in  the  scope  of  this 
research,  as  wings-level  flight  was  assumed. 

4.  All  objects  are  assumed  to  be  in  the  camera’s  field  of  view.  In  a  realistic 
scenario,  the  effects  of  wide  field  of  view  optics  or  multiple  optics  would  need 
to  be  considered. 

1.5  Thesis  Overview 

The  following  section  outlines  the  remainder  of  this  thesis.  Chapter  2  proceeds 
through  the  mathematical  basis  behind  the  research  and  discusses  alternative  cur¬ 
rent  methods  of  non-GPS  navigation  and  feature  extraction.  Chapter  3  presents  the 
methodology  of  how  the  research  was  conducted  and  the  specific  mathematical  mod¬ 
els  that  were  incorporated.  Chapter  4  shows  the  results  of  the  trajectory  estimation 
based  on  an  INS  aided  by  lunar  feature  angle  measurements.  Chapter  5  presents  the 
conclusion  and  suggestions  for  future  work  on  this  subject. 
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II.  Background 


2.1  Introduction 

The  following  section  discusses  the  theory  and  various  techniques  supporting  iner¬ 
tial  navigation  integrated  with  lunar  feature  angle  measurements.  Section  2.2  presents 
mathematical  notation  and  constructs.  Section  2.3  lays  out  the  frames  of  reference 
and  how  they  will  be  annotated.  Section  2.4  discusses  the  Kalman  Filter,  Extended 
Kalman  Filter,  and  error  state  estimation.  Section  2.5  covers  computer  vision  and 
camera  projection  theory.  Section  2.6  covers  the  Sun,  Moon,  and  Earth’s  orbital 
parameters  as  well  as  a  brief  description  of  the  lunar  surface  features.  Section  2.7 
presents  the  literature  review  covering  various  image-aided  navigation  techniques  and 
current  use  of  natural  passive  signals  to  gain  navigation  information.  Lastly,  Section 
2.8  summarizes  the  topics  discussed  in  this  Chapter. 

2.2  Arithmetic  Notation 

The  following  notation  will  be  used  for  the  rest  of  this  thesis: 

•  Upper  and  lower  case  italicized  letters  are  representative  of  scalar  variables. 

•  Bold  lower  case  letters  represent  arrays. 

•  Upper  case  rectified  letters  denote  a  cartesian  coordinate  axis. 

•  Bold  upper  case  letters  represent  matrices. 

—  I  is  an  identity  matrix  with  3  rows  and  3  columns. 

—  A  zero  0  in  bold  font  is  a  symbol  for  a  3  rows  by  3  columns  matrix  of  zeros. 

—  T  as  a  superscript  of  a  matrix  or  an  array  denotes  a  transpose,  i.e.  AT. 
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—  Direction  Cosine  Matrices  (DCM)  are  matrices  containing  information 
which  is  used  to  move  a  vector  from  one  frame  of  reference  to  another. 
DCMs  are  represented  in  the  matrix  notation  form  with  a  bold  upper  case 
letter  containing  a  superscript  and  a  subscript  variable  as  shown  in  Equa¬ 
tion  2.1  [28].  The  subscript  variable  n  is  the  current  frame  of  reference. 
The  superscript  variable  b  is  the  reference  frame  to  be  transitioned  to.  By 
transposing  a  DCM,  the  transition  from  one  frame  to  another  is  reversed 
as  in  Equation  2.2  where  the  transition  from  the  n-frame  to  the  b- frame  is 
reversed.  Section  2.3  provides  an  explanation  for  the  frames  of  reference 
used  in  this  research  effort. 

<7  (2.i) 

C£  =  (CS)T  (2.2) 

•  Skew  symmetric  matrices  are  matrices  populated  by  the  elements  of  a  single 
vector  or  array  [28].  An  array  variable  next  to  a  x  symbol  is  representative  of 
a  skew  symmetric  matrix.  Specifically,  in  this  research  effort,  a  skew  symmetric 
matrix  is  shown  as  a  3  rows  by  3  columns  matrix  formed  from  the  elements  of 
a  3  rows  by  1  column  vector  as  shown  in  Equations  2.3  and  2.4. 

7i 

7  =  y2  (2.3) 

73 

0  -73  72 

7X  =  73  0  — 7i  (2.4) 

-72  7i  0 
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2.3  Reference  Frames 


In  navigation,  distinct  frames  of  reference  are  needed  to  explain  an  aircraft’s  tra¬ 
jectory  with  respect  to  different  origin  points.  The  following  subsections  discuss  each 
frame  used  in  this  research  effort  as  described  in  [28]. 

Earth  Centered  Inertial  Frame  (i- frame). 

The  i-frame  is  a  cartesian  coordinate  reference  frame  with  origin  centered  at  the 
Earth’s  center  Op  The  three  axes  within  the  coordinate  system  are  non-rotating 
being  considered  fixed  to  the  stars  [28].  Figure  2.1  shows  the  three  axes  as  solid 
black  lines  with  axes  denoted  with  the  subscript  i.  The  Z-axis  is  pointed  through  the 
Earth’s  north  pole,  the  X-axis  is  pointed  through  the  Earth’s  equator  towards  the 
star  Aries,  and  the  Y-axis  points  perpendicular  to  the  X-axis  vector  along  the  Earth’s 
equator  plane  [28]. 

Earth  Centered  Earth  Fixed  Frame  (ECEF)  (e-frame). 

The  e-frame  is  a  cartesian  coordinate  reference  frame  with  origin  fixed  at  the 
Earth’s  center  Oe  as  shown  in  Figure  2.1.  The  e-frame  axes  are  denoted  by  a  subscript 
e.  The  Z-axis  lies  through  the  center  of  the  Earth  and  the  poles  shared  by  the  i-frame 
Z-axis.  The  X-axis  lies  perpendicular  to  the  Z-axis  through  the  center  of  the  Prime 
Meridian,  and  the  Y-axis  lies  at  longitude  90°E  pointing  through  the  Earth’s  equator 
[30].  The  X  and  Y  axes  rotate  around  the  Z-axis  at  the  Earth’s  rate  of  rotation  (0), 
which  is  7.292115  x  10“5^  (WGS-84)  [28], 

Earth  Fixed  Navigation  Frame  (n-frame). 

The  n-frame  is  a  three  dimensional  (3D)  cartesian  coordinate  reference  frame  that 
is  aligned  to  the  local  area  directions  of  north,  east,  and  down  (NED)  forming  a  plane 


tangent  to  the  Earth’s  surface  [30].  The  down  direction  corresponds  to  the  direction 
of  the  Earth’s  gravity  towards  the  center  of  the  Earth.  In  this  research  effort,  the 
orientation  of  the  n-frame  is  fixed  to  the  aircraft’s  trajectory  starting  position.  Figure 
2.1  illustrates  the  n-frame  axes  with  a  subscript  n.  The  n-frame  origin  is  at  On  along 
the  Earth’s  surface. 

Body  Frame  (b-frame). 

The  b-frame  follows  a  3D  cartesian  coordinate  system.  Figure  2.2  presents  an 
aircraft  body  with  the  origin  at  the  vehicle  body  center  Ob  [28].  The  body  frame 
axes  are  denoted  by  a  subscript  b.  The  X-axis  corresponds  to  the  roll  axis  and  points 
through  the  nose  of  the  vehicle  body.  The  Y-axis  corresponds  to  the  aircraft  pitch 
and  points  through  the  aircraft  wing  direction.  The  Z-axis  is  tied  to  the  aircraft  yaw 
and  points  vertically  through  the  center  of  the  aircraft. 

Sensor  Frame  (s-frame). 

The  s-frame  is  a  reference  frame  constructed  on  the  orientation  of  a  sensor  (CCD 
camera)  fixed  to  the  aircraft  body.  Figure  2.3  shows  the  sensor  with  the  origin  at  Os. 
For  this  research  effort,  the  sensor  will  be  fixed  to  the  b-frame  with  the  lens  and  Z-axis 
pointed  towards  the  sky  to  view  the  Moon.  The  X  and  Y  axes  point  perpendicularly 
to  each  other  along  a  plane  parallel  to  the  image  plane.  The  three  sensor  axes  are 
denoted  by  a  subscript  s. 
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Figure  2.1.  ECEF  and  Navigation  Frames.  Used  with  permission  [30]. 


Figure  2.2.  Vehicle  Body  Frame.  Used  with  permission  [30]. 
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Figure  2.3.  Sensor  Frame. 


2.4  Estimation 

Estimation  is  the  process  by  which  an  unknown  quantity  can  be  predicted  within  a 
margin  of  error  [20].  The  following  subsections  discuss  dynamic  systems,  the  Kalman 
Filter  (KF),  the  Extended  Kalman  Filter  (EKF),  and  error  state  estimation. 

Dynamic  Systems. 

Equation  2.5  represents  a  time  invariant  dynamic  system  in  the  form  of  a  linear 
differential  equation  characterized  by  continuous  time  [5]  [20].  The  dynamic  system 
models  state  behavior  as  it  is  propagated  over  time  [5].  A  state  is  a  phenomena  of 
interest,  e.g.  position,  velocity,  or  attitude  [5].  Fx(£)  in  Equation  2.5  represents  the 
states  mapped  to  the  dynamic  system  where  x(t)  is  the  state  vector  and  F  is  the 
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state  mapping  matrix.  Bu(£)  represents  the  system  inputs  mapped  onto  the  dynamic 
system  where  u (t)  is  the  system  input  vector  and  B  is  the  system  input  mapping 
matrix  [20].  Gw(t)  represents  the  system  noise  where  w(t)  is  the  noise  vector  and  G 
is  the  system  noise  mapping  matrix.  Equation  2.6  presents  Q (t)  as  the  input  noise 
covariance  matrix  where  the  Dirac  delta  function  is  represented  by  S(t)  [21]  [30]. 


x(f)  =  Fx(f)  +  Bu(t)  +  Gw  (t) 


(2.5) 


Q  (i)i(r)  =  E 


w(f)wT(f  +  r) 


(2.6) 


Equation  2.7  represents  the  state  uncertainty  where  P (t)  is  the  state  covariance  ma¬ 
trix.  The  state  covariance  propagation  equation  is  shown  in  Equation  2.8  [20]. 


Pxx(t) 


E 


x(f)xT(f) 


(2.7) 


Pxx(t)  =  FPxx(f)  +  Pxx(f)FT  +  GQ(f)GT  (2.8) 

Equations  2.9-2.12  represent  the  continuous  dynamic  system  from  Equations  2.5 
and  2.8  brought  into  discrete  form.  This  conversion  is  needed,  because  the  EKF 
algorithm  used  in  this  research  effort  operates  in  discrete  time  [5] .  For  the  rest  of  this 
thesis,  all  discrete  form  equations  hold  a  time  step  of  (A t)k  where  k  is  the  time  step 
index  number.  In  Equations  2.9  and  2.10,  is  the  state  transition  matrix  and  Q<j  is 
the  discrete  input  noise  covariance  matrix  respectively  [5]  [20].  Equation  2.11  shows 
the  discrete  propagation  of  the  state  mean  and  Equation  2.12  shows  the  discrete 
propagation  of  the  state  variance. 


<j>  =  g  F(At)fc 


(2.9) 
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Qd  =  E  wfcwjr 

(2.10) 

Xfc  =  4>Xfc_i  +  Wfc 

(2.11) 

k  =  ^Pfc-i^T  +  Qd 

(2.12) 

Kalman  Filter. 

The  KF  is  a  recursive  algorithm  that  performs  minimum  mean  square  error 
(MMSE)  estimation  of  state  dynamics  [5].  The  KF’s  algorithm  design  assumes  a 
Gaussian  distribution  for  the  state  dynamics  to  be  estimated.  This  assumption  al¬ 
lows  the  filter  to  completely  describe  a  state’s  dynamics  by  estimating  the  state’s 
mean  and  variance  or  uncertainty  [5] .  Three  components  comprise  the  KF :  the  prop¬ 
agation  phase  in  Equations  2.13  and  2.14,  the  measurement  model  in  Equation  2.15, 
and  the  update  phase  shown  in  Equations  2.18-2.20. 

Propagation  Phase. 

The  propagation  phase  is  comprised  of  propagating  the  state  mean  xy  and  variance 
P/,  estimates  through  discrete  time  steps.  When  no  measurement  update  is  present, 
the  propagation  equations  recursively  use  the  previous  state  mean  x^T  and  variance 
P;T  estimates  to  obtain  its  current  estimation  of  the  state  mean  and  variance  P^ 
as  shown  in  Equations  2.13  and  2.14  [5]. 


V  =  ***-i+w* 

(2.13) 

P*  =  4>Pti*T  +  Qd 

(2.14) 
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Measurements. 


KF  estimation  operates  through  the  introduction  of  sensor  measurements  z/,.  into 
the  filter  at  specified  discrete  times  [5].  The  measurements  aid  the  filter  by  provid¬ 
ing  information  about  the  states.  Equation  2.15  below  is  characteristic  of  a  linear 
measurement  model  [5] .  H/;.x/.  represents  the  relationship  between  the  sensor  and  the 
states  where  H*,  is  the  measurement  matrix.  The  level  of  information  to  which  the 
sensor  knows  about  the  states  is  dependent  upon  the  sensor  specifications  denoted  by 
the  state  observability  quantities  in  H/  [5] .  The  measurement  noise  v/,:  is  assumed  to 
be  white  Gaussian  noise.  Equation  2.16  is  the  measurement  noise  covariance  matrix 


Rfc- 


Zfc  H/jX/j  T  Vfc 


(2.15) 


Rfc 


E 


T 


(2.16) 


Update  Phase. 

The  final  component  of  the  KF  algorithm  is  the  update  phase.  During  the  update 
phase,  aiding  measurements  z/i:  are  input  into  the  algorithm  to  bring  about  an  optimal 
estimate  of  the  state  mean  and  variance.  The  criterion  held  to  determine  an  optimal 
estimate  within  the  KF  is  an  estimate  with  minimal  mean  square  error  P&  [5].  To 
achieve  this  aim,  an  estimator  based  upon  Bayes  rule  is  used  where  the  optimal 
estimate  x^  of  the  states  x*,  is  conditioned  on  the  measurement  z*,  [5] . 


xj  =  E(xt|zfc) 


(2.17) 


Based  on  the  above  Bayesian  estimator,  the  associated  updated  mean  and  covariance 
equations  are  shown  in  Equations  2.18-2.20  where  K&  is  the  Kalman  gain.  The 
Kalman  gain  is  a  blending  factor  that  works  to  blend  the  measurement  with  the  a 
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priori  state  estimate  to  bring  about  the  optimal  updated  estimate  [5]. 


=  S*  +  Kt(zi  -  Htxt ) 


(2-18) 


PJ  =  (I  -  K*H*)Pt- 

(2.19) 

(2.20) 

Extended  Kalman  Filter. 


A  modified  KF  algorithm  called  an  EKF  is  designed  to  process  non-linear  measure¬ 
ment  models  [5].  Equations  2.21-2.26  represent  the  EKF  algorithm  [5].  The  propaga¬ 
tion  phase  consists  of  Equations  2.21  and  2.22  which  are  the  estimated  state  mean  and 
variance  propagation  equations  respectively.  Equation  2.21  shows  the  states  propa¬ 
gated  by  a  non-linear  function  f(xfc)  [21].  The  EKF  covariance  propagation  equation 
in  Equation  2.22  holds  the  same  form  as  the  KF  covariance  propagation  equation  in 
Equation  2.14. 


(2.21) 


Pfc  =  ^-i^T  +  Qrf 


(2.22) 


The  EKF  update  phase  consists  of  Equations  2.23  and  2.24.  The  residual  holds  a 
non-linear  function  h(xfc)  in  Equation  2.23  within  its  calculation  [21].  The  Kalman 
gain  Kfc  in  both  equations  is  the  same  form  as  the  Kalman  gain  equation  in  the 
KF  algorithm  shown  in  Equation  2.20.  The  covariance  update  equation  in  Equation 
2.24  holds  the  same  form  as  the  covariance  update  equation  in  the  KF  algorithm  in 
Equation  2.19.  The  form  of  the  linear  measurement  matrix  H/r.  accommodating  the 
non-linear  measurement  function  will  be  presented  in  Equation  2.26. 


=xfc  +Kk(zk  -h(xj) 


(2.23) 
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Pfc  =  (I  -  K*H*)Pt- 


(2.24) 


The  measurement  model  for  an  EKF  is  likewise  non-linear  as  shown  in  Equation 
2.25.  The  linear  mapping  of  the  state  observabilities  onto  the  state  estimate  vector 
HfcXfc  within  the  KF  algorithm  is  replaced  by  passing  the  estimated  states  through  a 
non-linear  function  h(xfc)  which  holds  information  about  the  states  [5]  [21]. 


zk  =  h(xfc)  +  vfc  (2.25) 

To  process  the  non-linear  measurement  function  h(xfc)  through  the  EKF,  the  non¬ 
linearity  of  the  function  must  be  accounted  for  as  the  EKF  algorithm  requires  a 
linear  measurement  matrix  Hfc  when  calculating  the  Kalman  gain  Kfc  and  updated 
variance  estimate  as  shown  in  Equations  2.20  and  2.24  respectively  [21].  During 
the  update  phase,  a  method  is  used  to  linearly  approximate  around  the  non-linear 
measurement  function  h(xfc)  to  approximate  Tp.  A  linear  approximation  of  a  non¬ 
linear  function  can  be  found  through  the  use  of  a  first-order  Taylor  series  expansion 
known  as  a  Jacobian  which  consists  of  partial  differentiation  applied  to  the  non-linear 
measurement  function  as  seen  in  Equation  2.26  [21]. 

H  k  =  (xfc)  (2.26) 


Error  States. 

So  far,  modeling  and  estimation  has  been  discussed  with  regard  to  direct  states 
Xfc.  In  the  case  of  an  INS,  the  system  performance  is  characterized  by  the  amount  of 
trajectory  drift  error  over  time  [28].  Thus,  modeling  of  the  INS  error  states  Jxfc  are 
more  suitable  for  the  INS  design.  Also,  a  relationship  can  be  constructed  between  the 
errors  Jxfc,  errant  truth  xjns,*:,  and  true  states  Xfc  as  shown  in  Equation  2.27  [5] [21]. 
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Lastly,  the  error  states  can  be  modeled  in  the  same  way  as  the  direct  states  just  by 
modifying  the  state  vector  x/;.  to  incorporate  state  errors  as  shown  in  Equation  2.28 
[28]- 


Xfc  —  XiNS,fc  +  <5Xfc 

(2.27) 

5xk  =  $hxfc_i  +  wk 

(2.28) 

2.5  Computer  Vision 

This  section  provides  an  overview  of  basic  optical  sensor  configuration,  camera 
projection  theory,  angular  field  of  view,  and  angular  resolution  within  diffraction 
limited  systems. 

Camera  Projection  Theory. 

Figure  2.4  illustrates  a  pinhole  camera  model.  Light  reflected  or  projected  from 
objects  in  the  scene  of  the  outside  world  is  taken  in  by  the  aperture  (or  opening) 
shaped  as  a  double  convex  lens  into  a  dark  space  where  the  image  is  projected  [10]. 
Based  upon  the  behavior  of  light  rays  entering  the  aperture,  the  image  is  naturally 
an  inverted  projection  of  the  outside  scene  [10].  Depending  on  the  focal  length  /, 
the  projected  arrow  image  can  be  manipulated.  A  large  focal  length  will  create  a 
magnified  image  of  the  scene  with  a  small  field  of  view.  A  small  focal  length  will 
cause  the  scene  to  be  viewed  from  a  wide  field  of  view.  Equation  2.29  shows  the 
relationship  between  the  scene,  image,  and  focal  length  where  Si  is  the  distance  from 
the  scene  to  the  aperture  and  S2  is  the  distance  from  the  lens  to  the  projected  scene 
[10]- 
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SCENE 


>k - s2 - >j 


IMAGE 


Figure  2.4.  Pinhole  Camera  Model.  Used  with  permission  [30]. 


1  1  _  1 

Si  s2  7 


(2.29) 


ILLUMINATION 


Figure  2.5.  Optical  Sensor  Model.  Used  with  permission  [30]. 

Figure  2.5  follows  the  pinhole  camera  concepts  presented  in  Figure  2.4,  but  builds 
upon  the  model’s  analog  aspects  with  an  analog-to-digital  interface  converting  raw 
images  of  the  outside  world  scene  into  digital  information  a  computer  can  use  to 
reconstruct  the  image.  Within  Figure  2.6,  the  image  projection  from  Figure  2.4 
is  placed  in  the  context  of  the  3D  outside  scene  projected  onto  a  2D  image  plane. 
In  Figure  2.6,  sc  is  the  object  location  in  the  outside  scene  and  spro’  is  the  object 
location  on  the  2D  projection  [30].  The  2D  image  axes  are  represented  by  a  2D 
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cartesian  coordinate  system  with  xproj  and  ypr0j  as  the  horizontal  and  vertical  axes 
respectively.  The  relationship  between  the  object  location  in  the  3D  scene  sc  and 
the  object  projection  on  to  the  image  spro>  is  shown  in  Equation  2.30  where  zc  is  the 
direction  from  the  camera  to  the  3D  scene. 


Figure  2.6.  Camera  Projection  Model.  Used  with  permission  [30]. 


spr°j  =  JL  sc 


(2.30) 


The  object  projection  location  needs  to  be  determined  with  respect  to  pixel  loca¬ 
tions  in  the  digital  image  as  shown  in  Equation  2.31  [30].  M  and  N  serve  as  the  digital 
array  horizontal  and  vertical  axes  respectively  of  width  W  and  height  H.  Figure  2.7 
shows  the  2D  digital  image  array  coordinate  system  relating  the  projected  object 
location  spro>  coordinates  onto  the  digital  image  array  coordinate  system  [30]. 


M 

H 

0 

0 

spr°j  + 

Af+1 

2 

0 

N 

W 

0 

N+l 

2 

(2.31) 
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Figure  2.7.  Digital  Image  Coordinate  System.  Used  with  permission  [30]. 


Angular  Field  of  View. 

The  angular  field  of  view  is  the  angle  of  the  world  scene  observed  by  a  camera’s 
lens  [10].  For  a  camera  with  a  rectangular  image  array,  the  resulting  angular  field  of 
view  is  described  by  two  angles:  a  width  angle  aw  and  a  height  angle  [1].  Figure 
2.8  presents  a  side  view  of  a  camera  model  showing  only  the  height  field  of  view 
parameters.  Equation  2.32  shows  the  height  angle  field  of  view  equation  with  respect 
to  the  image  array  height  /if  and  focal  length  [1]  [10] . 

ah  =  2tan_1  (J^j  (2-32) 

Equation  2.33  shows  the  determination  of  the  height  of  the  world  scene  hscene  observed 
by  the  camera  where  dscene  is  the  distance  from  the  lens  to  the  world  scene. 


h 


scene 


2d 


scene 


tan 


(2.33) 
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To  determine  the  angular  field  of  view  for  width  crw,  Equation  2.32  can  be  used  and 
the  image  array  height  replaced  by  the  image  array  width  Wf. 


World  Scene 


Figure  2.8.  Angular  field  of  view  based  on  image  array  height  (Side  View). 


Diffraction  Limited  Spot  Size. 

This  subsection  discusses  the  resolution  of  a  point  source  within  an  image  con¬ 
strained  by  lens  diffraction.  In  ideal  circumstances,  a  lens  would  perfectly  transfer  a 
point  source  in  the  world  scene  to  its  detected  location  on  an  image  array  as  shown 
in  Figure  2.9  [10].  A£  denotes  the  center  to  center  pixel  separation  or  pixel  size. 

Point  Source 


Figure  2.9.  Ideal  Resolution  of  Point  Source  Along  Image  Array. 
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Unfortunately,  due  to  diffraction  caused  by  a  camera’s  finite  aperture  (neglecting 
atmospheric  effects),  the  light  of  the  point  source  is  spread  over  the  detector  [10].  The 
diffraction  limitation  of  the  lens  causes  angular  variance  A 6q  in  the  location  of  the 
point  source  projection  within  the  pixels  [10].  The  resulting  projection  of  the  point 
source  on  the  image  results  in  an  Airy  disk  or  diffraction  limited  spot,  meaning  that 
the  point  source  projection  can  be  located  somewhere  within  the  first  ring  of  the  Airy 
disk  radius  qi  [10].  Figure  2.10  showcases  the  angular  resolution  of  a  point  source  in 
the  world  scene  upon  pixel  squares  within  an  image  distorted  by  diffraction. 


World  Scene 


Point  Source 


Figure  2.10.  Angular  Resolution  of  Point  Source  Along  the  Image  Array. 


The  conversion  between  qi  and  A 0q  is  shown  in  Equation  2.34  [10].  Equations  2.35 
and  2.36  describe  the  relationships  between  the  angle  of  resolution  A 9q  and  Airy  disk 
radius  q\  to  the  focal  length  /,  light  wavelength  Aq,  and  lens  diameter  D  respectively 
[10]. 

Adq  ~  sinA^q  =  y  (2.34) 

^  (2-35) 
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1.22/Aq 
ql  =  ^~ 


(2.36) 


Depending  on  the  pixel  size  A£,  a  mismatch  in  sampling  could  occur  [10].  Figure 
2.11  presents  three  cases  of  sampling.  Case  A  shows  an  over  sampled  case  where  the 
Airy  disk  size  is  spread  over  multiple  pixels.  Case  B  presents  the  critically  sampled 
case  where  the  sampling  size  is  equal  to  the  Airy  disk  diameter  cordoning  the  point 
source  projection  to  one  pixel.  Case  C  is  an  under  sampled  case  where  the  pixel 
size  is  much  larger  than  the  point  source  projection  creating  ambiguity  to  where  the 
projection  is  located  within  the  pixel  itself.  To  obtain  sub-pixel  resolution  of  the  Airy 
disk,  it  is  preferred  to  have  the  Airy  disk  spread  over  many  pixels  [10]  [18]  . 


Figure  2.11.  Angular  resolution  of  diffraction  limited  spots  with  respect  to  sample  size 
along  the  image  plane. 

2.6  Overview  of  Moon,  Sun,  and  Earth  Orbital  Parameters 

The  following  section  provides  a  high  level  overview  of  the  celestial  orientation 
between  the  Sun,  Moon,  and  Earth.  Table  2.1  presents  the  key  parameters  which 
summarize  the  relationship  between  the  three  bodies. 

Figures  2.12-2.14  illustrate  the  Earth  holding  an  axial  tilt  of  23.44°  with  respect  to 
its  orbit  plane  around  the  Sun  [31].  The  Earth’s  sidereal  orbit  period  around  the  Sun 
is  365.256  days  traveling  in  a  counter  clockwise  direction  [31].  The  Earth  performs  a 
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Table  2.1.  Sun  [32],  Moon  [33],  and  Earth  [31]  Parameters. 


Celestial  Body 

Radius  (km) 

Distance  to 

Earth  (km) 

Inclination  to 
Equator  (°) 

Inclination  to 
Ecliptic  (°) 

Obliquity  to 
Orbit  (°) 

Sun 

696,000 

149,600,000 

N/A 

N/A 

N/A 

Moon 

1,738.1 

378,000 

18.28-28.58 

5.145 

6.7 

Earth 

6,378.1 

N/A 

N/A 

0 

23.44 

revolution  around  its  own  axis  every  23.93  hours  [31].  The  ecliptic  is  the  plane  of  the 
Earth’s  orbit  around  the  Sun  [31]. 


Figure  2.12.  View  of  Orientation  of  Earth,  Moon,  and  Sun  Orbits  [31]. 


The  Moon’s  revolution  period  around  the  Earth  is  27.32  days  [33].  Figure  2.14 
presents  the  Moon’s  orbit  inclination  with  respect  to  the  Earth’s  equator  ranging  from 
18.28°  to  28.58°  [33].  During  the  Moon’s  orbit  around  the  Earth,  the  Moon’s  axial 
tilt  varies  in  librations  along  its  latitude  and  longitude.  Latitude  librations  occur  at 
6.7°,  and  longitude  librations  occur  at  7.7°  [11]. 
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Figure  2.13.  Top  View  of  Orientation  of  Earth  [31],  Moon  [33],  and  Sun  [32]  Orbits. 


Figure  2.14.  Side  View  of  Orientation  of  Earth  [31],  Moon  [33],  and  Sun  [32]  Orbits. 
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Figure  2.15  presents  a  high  resolution  image  of  the  Moon’s  near  side  captured  by 
the  Lunar  Reconnaissance  Orbiter  wide  angle  camera  (LROC  WAC)  in  partnership 
between  NASA  and  Arizona  State  University.  The  lunar  near  side  surface  contains 
unique  characteristics  defined  by  smooth  basins  called  seas  or  mare,  craters,  channels, 
mountains,  ejecta  rays,  and  ridges  [11].  The  lunar  mare  are  the  dark  smooth  areas 
along  the  Moon’s  surface  in  Figure  2.15  [11].  Ejecta  rays  are  the  displaced  materials 
of  impact  craters  that  have  fallen  into  place  radially  around  the  impact  crater  [11]. 
Ejecta  rays  can  be  seen  in  Figure  2.15  as  the  lighter  straight  lines  oriented  radially 
around  the  crater  named  Mount  Copernicus  in  the  middle  to  left  portion  of  the  image 
lying  surrounded  by  the  dark  mare  (south  of  Mare  Ibrium)  [11], 


Figure  2.15.  Lunar  Near  Side.  Used  with  Permission  [26]. 

2.7  Literature  Review 

This  section  discusses  current  non-GPS  navigation  sensors  found  in  open  literature 
which  extract  navigation  information  from  natural  phenomena.  The  first  subsection 
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presents  attitude  determination  by  computer  vision  techniques  based  upon  lunar  and 
planetary  features.  The  second  subsection  discusses  other  various  attitude  determina¬ 
tion  techniques  derived  from  knowledge  of  the  Sun’s  position,  Earth’s  magnetic  field, 
and  Earth’s  horizon,  etc.  Attitude  refers  to  a  vehicle’s  orientation  in  the  b-frame 
represented  by  the  Euler  angles:  roll  0,  pitch  9,  and  yaw  ip  [28].  Aiding  attitude 
updates  have  been  shown  to  improve  the  uncertainty  of  KF  estimated  INS  position 
error  states  [2], 

Attitude  Determination  by  Lunar  or  Planetary  Features. 

Three  processes  in  lunar  or  planetary  image-aided  navigation  methods  are:  image 
collection,  navigation  data  extraction,  and  attitude  determination  [7].  Planetary  and 
lunar  images  have  been  identified  as  an  attitude  information  source  with  respect  to 
a  vehicle’s  body  frame  [6]  [24] . 

In  space,  a  vehicle’s  body  frame  attitude  could  be  found  with  respect  to  a  reference 
frame  based  on  a  plane  consisting  of  the  vehicle  and  two  celestial  bodies  denoted  by 
the  p-frame  in  Figure  2.16.  The  p- frame  origin  Op  and  the  b-frame  origin  of  the  vehicle 
Ob  share  the  same  location.  The  pitch  9  and  yaw  ip  angles  of  the  b-frame  with  respect 
to  the  p-frame  can  be  estimated  based  on  the  plane  established  by  the  three  bodies 
as  shown  in  Figure  2.16  [8].  The  plane  in  Figure  2.16  consists  of  the  Moon,  Sun, 
and  a  vehicle  with  the  plane’s  origin  at  Op.  The  roll  angle  <p  of  the  vehicle  body 
with  respect  to  the  p-frame  is  estimated  based  upon  the  measured  rotation  angle  9 
of  the  vehicle  body  in  relation  to  the  Moon.  Mortari  and  Park  propose  using  a  stored 
catalog  Moon  image  to  compare  with  the  observed  Moon  image.  This  would  allow  a 
translation  and  rotation  calculation  between  the  images  to  be  performed  to  determine 
0  [24],  Enright’s  method  to  determine  t)  differed  due  to  his  use  of  over-exposed  star 
tracker  Moon  images  instead  [8]. 
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Figure  2.16.  Attitude  Determination  by  Celestial  Bodies  [8]. 

Key  Points  Feature  Detection. 

Key  points  feature  detection  and  extraction  is  characterized  by  the  detection  and 
tracking  of  features  of  interest.  Absolute  attitude  sensors  such  as  second  generation 
star  trackers  use  the  method  to  detect  and  track  celestial  bodies,  i.e.  stars  [18].  The 
method  can  also  be  used  on  planetary  features  [15].  The  process  works  by  having 
key  points  information  stored  within  a  processing  software  database  or  catalog.  The 
database  holds  information  about  each  key  point’s  unique  characteristics  and  can  be 
compared  with  an  observed  image  of  the  same  terrain  or  star  pattern  to  gain  rotation 
and  translation  information  in  a  vehicle’s  body  frame  [7] . 

For  the  star  tracker,  a  star’s  magnitude  and  brightness  serve  as  unique  qualities 
and  are  stored  in  a  star  catalog  [16].  The  star  catalog  holds  information  about 
each  stars  arrangement  in  relation  to  the  nearest  stars.  A  set  of  three  stars  in  close 
proximity  is  called  a  Triangle  Feature  (TF)  [16].  This  can  be  seen  in  Figure  2.17.  TFs 
are  the  basic  pattern  from  which  meaningful  attitude  information  can  be  extracted 


28 


within  a  star  tracker  [16].  The  unique  characteristics  of  a  TF  are  listed  in  Equation 
2.37  derived  from  Figure  2.17  where  Star  C  is  the  center  star,  Star  A  is  the  first 
nearest  star  to  Star  C,  and  Star  B  is  the  second  nearest  neighbor  to  Star  C.  Within 
[16],  the  star  catalog  stores  a  vast  number  of  TF  combinations.  Liebe  in  [16]  performs 
a  quantization  on  the  TF  combinations  to  conserve  memory.  Under  quantization,  the 
stored  TF  reference  listed  in  Figure  2.17  would  read  [55,  81,  32].  A  star  averages  180 
star  catalog  references  in  relation  to  the  TFs  its  contained  in  [16]. 

Distance  to  Star  A,  Distance  to  Star  B,  Angle  between  Stars  A  and  B 


(2.37) 


Figure  2.17.  Star  Tracker  Triangle  Feature  [16]. 

In  many  cases,  it  is  not  enough  to  locate  key  points  based  at  given  pixel  locations. 
Subpixel  resolution  quality  is  needed  to  determine  exactly  where  the  key  point  is 
located  within  an  image.  This  is  called  hyperacuity  or  subpixel  precision  [18].  To 
achieve  this  goal,  the  process  of  centroiding  can  be  used.  At  the  pixel  level,  a  key 
point  can  appear  to  be  uniform  in  intensity.  Centroiding  consists  of  slightly  blurring 
an  image  to  obtain  a  Gaussian  form  of  pixel  intensities  around  a  key  point  [18] .  The 
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tighter  mean  extracted  from  the  key  point  intensity  distribution  provides  a  more 
accurate  subpixel  level  data  point  supporting  an  improved  attitude  determination 
[18]- 

Van  Bezooijen’s  key  point  extraction  method  doesn’t  solely  rely  on  TFs  for  cor¬ 
respondence  between  the  observed  key  points  and  its  key  point  catalog  [7].  Van 
Bezooijen’s  algorithm  works  from  the  premise  that  an  increase  of  stars  in  view  cre¬ 
ates  a  more  unique  star  pattern  to  improve  attitude  determination.  Based  upon  Van 
Bezooijen’s  algorithm,  researchers  in  [7]  are  working  on  the  statistics  and  viability  of 
an  N-star  pattern  recognition  algorithm  for  attitude  determination. 

Key  points  extraction  and  detection  methods  can  be  used  in  extracting  planetary 
features  for  navigational  purposes  [15]  [17].  Liebe’s  criteria  in  [15]  for  a  robust  feature 
extraction  algorithm  upon  planet  terrain  is  for  the  algorithm  to  operate  proficiently 
withstanding  changes  in  vantage  points,  occlusion,  lighting,  and  rotation.  He  reaches 
the  conclusion  that  this  aim  can  be  achieved  through  the  use  of  token-based  closed 
contour  tracking  very  similar  to  key  points  extraction  and  detection  used  in  star  track¬ 
ers.  According  to  Liebe,  the  features  must  hold  a  distinguishable  quality  that  allows 
them  to  be  recognized  from  a  wide  range  of  vantage  points  and  lighting  conditions 
[15].  From  images  of  Jupiter,  he  was  able  to  extract  features  to  be  used  for  a  closed 
contour  feature  constellation  while  rotation  and  occlusion  separately  occurred  [17]. 

Hough  Transform. 

Another  method  of  feature  detection  and  extraction  that  has  been  explored  for 
attitude  determination  purposes  is  the  Hough  Transform.  The  Hough  Transform  is 
a  robust  algorithm  used  to  detect  specified  shapes  within  images,  particularly,  lines 
and  circles  [6]  [14] .  In  [6] ,  a  circular  crater  on  Saturn’s  moon,  Mimas,  was  successfully 
detected  using  the  Hough  Transform.  Both  [6]  and  [14]  voice  that  although  the 
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transform’s  performance  is  high  at  detecting  shapes  within  an  image,  an  increase  in 
the  parameters  and  complexity  of  the  detectable  image  shapes  would  greatly  increase 
the  computational  memory  needed.  In  [14],  the  suggestion  is  made  that  extensions  to 
the  Hough  Transform  exist  which  reduce  the  inefficient  algorithm  needs.  However,  the 
stripping  away  of  the  algorithm  leads  to  instability  in  the  algorithm  when  processing 
shapes  slightly  different  than  the  specified  template  [14]. 

Gray  Level  Cooccurance  (GLC)  Matrices. 

In  [29],  a  feature  detection  method  utilizing  GLC  matrices  was  implemented  to 
extract  planetary  features,  e.g.  mountains,  volcanoes,  and  valleys.  The  research  in 
[29]  was  conducted  in  the  context  of  providing  extracted  features  for  planetary  ter¬ 
rain  tracking  and  navigation  purposes.  A  GLC  matrix  scans  a  defined  image  region 
and  estimates  the  second-order  probability  for  a  specified  gray  level  change  between 
adjoined  pixels  [29].  Each  GLC  is  tuned  to  extract  a  type  of  image  feature.  Equation 
2.38  represents  a  type  of  GLC  tested  in  [29]  where  the  gray  levels  are  represented 
by  Lg,  6g  is  the  separation  between  any  two  pixels  within  the  local  region  Tg,  and  ig 
and  jg  are  the  gray  levels.  The  Inertia  GLC  matrix  below  suited  identifying  random 
gray-leveled  areas,  and  the  Entropy  GLC  matrix  was  found  to  be  the  most  robust 
GLC  matrix  [29] .  The  Cluster  Shade  GLC  matrix  was  found  to  favor  feature  extrac¬ 
tion  within  high  contrasting  gray  level  regions  which  are  characterized  by  craters. 
Although,  not  holding  high  performance  as  the  specialized  GLCs  in  their  respec¬ 
tive  environments,  it  was  consistent  in  its  performance  over  various  local  regions  and 
features. 

Inertia: 

Lg-lLg-l 

/(hg,Tg)=  y  ]  y  ]  (ig  —  jg)  s(ig,jg,  Sg,Tg)  (2.38) 

*g=°  ig=° 
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Following  the  GLC  matrix  tests  conducted  on  images  of  moons,  Triton  and  Mi¬ 
randa,  Udomkesmalee  and  Antonio  concluded  that  although  GLC  matrices  provided 
distinct  feature  detection  results,  the  quality  and  the  method  fell  short.  This  was 
due  in  part  to  the  difficulty  of  sorting  out  shadow  areas  causing  mis-characterizations 
of  planetary  features  and  the  complexity  surrounding  the  recognition  of  numerous 
feature  types  [29]. 


Template  Matching. 

Template  Matching  is  another  method  used  in  feature  detection  of  planetary  ter¬ 
rain  and  attitude  determination  as  read  in  [4]  and  [6].  In  [24],  a  method  of  tem¬ 
plate  matching  is  performed  on  Moon  images.  The  algorithm  called  the  Moon-Sun 
Attitude  Sensor  works  by  matching  an  observed  image  of  the  Moon  to  a  stored  ref¬ 
erence  image.  The  resulting  captured  image  is  then  translated  and  rotated  until  it 
matches  the  reference  Moon  image.  Equation  2.39  characterizes  the  match  of  the  two 
Moon  images  through  the  quantity  r  which  represents  the  correlation  quantity  [24], 
Once  r  equals  0,  then  a  perfect  match  has  been  reached.  Once  the  perfect  match  is 
achieved,  the  phase  angle  d,  which  the  template  was  rotated  to  match  the  observed 
image,  can  provide  information  to  aid  in  attitude  determination.  This  was  briefly 
discussed  in  the  beginning  of  this  section  and  shown  as  the  end  result  of  Mortari  and 
Park’s  attitude  algorithm  in  Figure  2.18.  In  Equation  2.39,  the  image  correlation 
equation  is  presented  where  Im  and  Jm  represent  the  respective  Moon  images  to  be 
matched.  7m  and  Jm  are  stated  to  be  the  mean  value  of  each  image,  c  and  d  are  the 
pixel  index  numbers. 


r  = 


Sc,d(7 n,(c,d)  Jm)  (Jm,(c,d)  Jm) 

(7m,(c,d)  —  Jm)  ('7n,(c,d)  —  Jm) 


(2.39) 
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A  similar  equation  to  Equation  2.39  was  found  in  [6]  showing  the  commonality  of  the 
correlation  equation  within  the  template  matching  method  of  feature  detection.  To 
increase  the  efficiency  of  correlation,  [6]  suggests  that  correlation  can  be  carried  out 
in  the  Frequency  Domain  with  Fast  Fourier  Transform. 


Figure  2.18.  Mortari  and  Park’s  Image  Processing  Outline.  Used  with  permission  [24]. 

Figure  2.18  is  a  feature  detection  algorithm  example  involving  template  matching. 
From  this  figure,  the  various  stages  of  image  processing  involved  in  template  matching 
feature  detection  can  be  viewed.  Mortari  and  Park  first  perform  median  filtering  to 
remove  noise  [24],  From  here,  a  binary  image  can  be  created  from  the  observed 
image  and  edge  detection  conducted.  As  unique  to  [24],  if  the  Moon  is  full  and 
symmetric,  the  sea  of  Crisium  will  be  selected  as  to  perform  edge  detection  upon.  At 
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this  point,  the  captured  edge  is  translated  and  rotated  to  gain  the  match  and  phase 
angle  i)  as  mentioned  earlier  in  this  section.  Lastly,  Mortari  admits  that  the  Moon- 
Sun  Attitude  Sensor  performance  would  not  match  the  performance  of  modern  star 
trackers;  however,  the  sensor  could  be  viable  in  budget  constrained  environments 
[23].  Also,  a  Moon-Sun  Attitude  Sensor  aided  INS  system  could  have  the  potential 
to  match  a  star  tracker’s  performance. 

Natural  Passive  Signals  Sensor  Methods. 

The  following  subsection  covers  different  sensor  types  which  derive  position  and 
attitude  information  from  natural  phenomena.  Each  listed  sensor  does  not  provide  a 
complete  navigation  solution  on  its  own,  but  provides  separate  angle  or  range  mea¬ 
surements  when  coupled  together  can  be  used  in  trajectory  estimation  or  to  gain  a 
position  fix.  Sensor  combination  addresses  this  issue  by  combining  the  measurements 
of  disparate  sensors  to  form  a  complete  navigation  solution  [13]. 

Sun  Sensor. 

The  Sun  holds  a  predictable  ephemeris  track  across  the  Earth  each  day.  From 
this  constant  motion,  angle  measurements  can  be  derived  from  different  Earth  loca¬ 
tions  with  respect  to  the  Sun  and  the  horizon  to  gain  an  elevation  angle  measurement 
which  holds  latitude  position  information  [9] .  Modern  sun  sensors  also  follow  a  similar 
elevation  angle  angle  concept,  however,  using  the  Sun’s  angles  for  attitude  determi¬ 
nation  instead  [18]  [19].  Liebe  has  continued  research  into  enhancing  the  ancient  sun 
sensor  methodology  to  develop  a  micro  sun  sensor  providing  solar  angle  information 
based  around  MEMS  technology  and  a  centroiding  technique  [19] .  Liebe’s  micro  sun 
sensor  consists  of  a  miniature  enclosed  box  with  an  aperture  at  the  top  to  create  a 
well  where  the  Sun’s  rays  can  enter.  Based  upon  the  direction  that  the  Sun’s  rays 
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strike  the  wafers  aligned  along  the  well  floor,  an  angle  of  the  Sun’s  rays  can  be  ob¬ 
tained  which  leads  to  a  unique  attitude  solution.  Liebe  in  [18]  lists  sun  sensors  with 
an  accuracy  of  1  arcminute  for  attitude  determination. 

Magnetometer . 

Magnetometers  are  sensors  which  employ  the  Earth’s  magnetic  field  to  gain  mag¬ 
netic  measurements  corresponding  to  a  unique  attitude  determination.  The  readings 
consist  of  the  size  and  orientation  of  the  magnetic  field  lines  [18].  The  reported 
accuracy  for  magnetometers  is  around  1  arcminute  [18]. 

Horizon  Sensor. 

Utilizing  measurements  based  upon  the  limb  of  the  Earth,  navigation  information 
can  be  derived  to  gain  a  unique  attitude  determination  [18] .  The  limb  measurements 
are  typically  captured  through  an  infrared  camera  [18].  Liebe  lists  the  accuracy  of 
horizon  sensors  at  5  arcminutes  which  is  the  least  accurate  as  compared  to  the  other 
sensors  listed  in  [18]. 

Star  Tracker. 

Referenced  in  the  previous  subsection,  star  trackers  are  optical  sensors  which 
capture  constellation  information  to  gain  an  attitude  determination.  In  this  research 
effort,  second  generation  star  trackers  are  referred  to.  Second  generation  star  trackers 
consist  of  a  charge  coupled  device  (CCD)  camera  and  an  associated  software  package 
which  holds  a  star  database  and  star  catalog  [16]  [18].  The  captured  constellation 
images  are  compared  to  a  stored  catalog  image  of  the  observed  constellation,  and  a 
unique  set  of  angles  can  be  computed  which  provide  a  navigation  solution  [18].  Star 
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trackers  hold  an  accuracy  of  1  arcsecond  as  shown  in  [18] .  This  is  the  highest  accuracy 
in  comparison  to  the  other  sensors  listed. 

2.8  Summary 

Chapter  2  provided  the  necessary  theoretical  background  to  support  the  methodol¬ 
ogy  approach  discussed  in  Chapter  3.  Dynamics  models,  computer  vision  techniques, 
and  camera  projection  theory  were  discussed  as  well  as  KF  and  EKF  estimation  meth¬ 
ods.  An  overview  of  the  lunar  orbital  parameters  and  the  lunar  surface  was  covered. 
Various  open  literature  methods  of  non-GPS  navigation  sensors  were  presented  as 
well. 
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III.  Methodology 


3.1  Introduction 

The  methodology  chapter  discusses  the  models  and  methods  that  were  put  into 
practice  to  perform  drift  correction  on  INS  errors  aided  by  lunar  feature  angle  mea¬ 
surements.  Section  3.2  discusses  the  overall  system  model.  Section  3.3  describes  the 
truth  model  used  to  simulate  the  flight  trajectory.  Section  3.4  presents  the  error  state 
motion  model  that  was  used  to  simulate  INS  errors.  Section  3.5  shows  the  validation 
of  the  error  state  motion  model.  Section  3.6  covers  the  measurement  model  used  to 
correct  the  INS  errors.  Section  3.7  presents  the  lunar  feature  model.  Section  3.8 
discusses  the  EKF  implementation  process.  Section  3.9  provides  the  summary  of  the 
section. 

3.2  System  Overview 

To  gain  an  understanding  of  the  system  as  a  whole,  the  following  section  discusses 
the  error  state  model,  the  error  state  estimation  model,  and  presents  the  system  block 
diagram. 

The  INS  error  model  in  Equation  3.1  shows  the  relationship  between  the  INS  error 
states  (jxfc,  errant  truth  xINg,fc,  and  true  states  x*.  [5].  The  true  states  x*.  hold  direct 
trajectory  information  and  the  trajectory  model  is  presented  in  Section  3.3.  The  INS 
error  states  foe*,  are  based  on  the  Pinson  15  error  state  dynamics  model  described  in 
Section  3.4.  Equation  3.2  shows  how  the  subtraction  of  the  INS  errors  from  the  true 
states  is  used  to  create  the  simulated  INS  trajectory  xiNs,fe-  With  the  error  states 
modeled  on  INS  error  noise  strengths  (see  Table  3.1),  the  corrupted  truth  can  act  as 
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a  realistic  INS  modeled  trajectory. 


Xfc  —  XiNS.fc  + 


(3.1) 


xiNS,fc  —  Xfc  —  (3.2) 

With  the  INS  trajectory  and  INS  errors  simulated,  estimation  of  the  true  states  x^, 
can  be  performed.  This  research  effort  utilizes  an  indirect  truth  estimation  approach 
called  error  state  estimation  whereby  the  INS  error  states  Sxk  are  estimated  as  shown 
in  discrete  form  cix/,.  in  the  error  estimation  propagation  equation  below  [5] . 

=  $5x+_1  +  wk  (3.3) 


Once  the  INS  error  states  are  estimated,  Equation  3.4  shows  the  quantities  are  added 
to  the  errant  INS  trajectory  xiNs,fc  to  bring  about  an  estimation  of  the  true  states 

Xjfe. 

xfc  =  xiNS;fe  +  hxfc  (3.4) 

Figure  3.1  presents  the  system  block  diagram  containing  the  INS  trajectory,  lunar 
vision  measurements,  and  the  EKF  working  together  to  estimate  an  aircraft’s  true 
trajectory  based  upon  the  error  state  estimation  model.  The  system  follows  an  indi¬ 
rect  feed  forward  format  [20]. 


Figure  3.1.  Lunar  Vision  Aided  Navigation  with  INS  Block  Diagram. 


3.3  Truth  Model 


The  truth  model  of  aircraft  motion  is  simulated  by  a  circular  trajectory  with 
radius  20km  along  the  north  and  east  plane  in  the  n-frame  as  shown  in  Equations 
3. 5-3. 8  with  distance  in  meters,  /i  is  the  circle  central  angle.  The  truth  along  the 
north  and  east  position  states  was  kept  unchanged  to  act  a  baseline  throughout  the 
various  scenarios  that  were  conducted,  Pn,&  and  vN  ^  are  the  position  and  velocity  of 
the  aircraft  along  the  north  axis  respectively. 


VN,fc 


Pn,o  +  20000sin(/r) 

(3.5) 

_  PN,fc  —  PN,/c-l 

(A  t)k 

(3.6) 

Pe,o  +  20000cos(/r) 

(3.7) 

_  PE,fc  —  PE,fc-l 

(A  t)k 

(3.8) 

VE,fc 

The  down  axis  motion  is  generated  by  a  First  Order  Gauss  Markov  (FOGM)  accel¬ 
eration  model  which  represents  aircraft  motion  in  a  realistic  fashion  [5].  Figure  3.2  is 
the  FOGM  acceleration  model  block  diagram  representing  the  down  axis  of  motion 
over  continuous  time,  r  is  the  system’s  correlation  time  constant  [20]. 


WD(t) 


aD(t)  vD(t)  pD(t) 


Figure  3.2.  n-frame  Down  axis  FOGM  acceleration  model  Block  Diagram  [20]. 
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The  set  of  first  order  linear  differential  equations  describing  the  down  axis  motion 
are  shown  in  Equations  3.9  through  3.11  where  each  differential  equation  is  defined 
by  the  motion  states  of  acceleration  a,  velocity  v,  or  position  p  [20].  w (t)  is  white 
Gaussian  noise  on  the  acceleration  state. 


a  d(£) 


y~~j  aD(f)  +  wD(f) 

(3.9) 

vd  (t)  =  a  D(f) 

(3.10) 

Pd(£)  =  vD(t) 

(3.11) 

The  first  order  differential  equation  form  of  the  down  axis  truth  model  is  shown  in 
Equation  3.12  with  incorporated  random  noise  [20].  However,  due  to  research  scope, 
the  control  inputs  u (t)  are  not  being  evaluated  in  the  state  propagation.  Therefore, 
the  control  inputs  are  set  to  zero  in  Equation  3.13. 


x(f)  =  Fx(/)  +  Bu(/)  +  Gw  (t) 


(3.12) 


x(t)  =  Fx(t)  +  Gw  (t) 


(3.13) 


Equation  3.14  represents  the  down  axis  state  space  truth  model. 


Pd 

0 

1 

0 

Pd 

0 

vd 

= 

0 

0 

1 

vd 

+ 

0 

aD 

0 

0 

_  1 

T 

aD 

wD 

(3.14) 


Equation  3.15  is  the  down  axis  input  noise  covariance  matrix  Q (f)  representing  the 
input  noise  for  the  down  axis  state  space  truth  model.  The  input  noise  is  white 
Gaussian  noise  with  zero  mean  [20].  The  relationship  between  the  input  noise  and 
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the  process  noise  can  be  seen  in  Equation  3.16  where  <td  is  the  process  noise  standard 
deviation  [20]. 

Q (t)  =  E  wD(t)w ■£(£)  =  E  w2  (3.15) 

Q  (t)  =  —  (3.16) 

T 

3.4  Error  State  Motion  Model 

To  simulate  realistic  INS  errant  behavior,  a  dynamics  model  following  the  con¬ 
tinuous  time  first  order  linear  differential  equation  form  with  incorporated  random 
noise  input  was  implemented  as  shown  in  Equation  3.17  [28].  The  control  inputs 
u(t)  are  not  being  evaluated  in  the  error  state  propagation  and  are  set  to  zero  in 
Equation  3.18.  The  error  states  are  seen  in  Equation  3.19  which  is  the  error  state 
vector  5x.  The  five  error  state  categories  are  position  error  5p,  velocity  error  5v,  tilt 
-0,  accelerometer  bias  5a,  and  gyroscope  bias  5b.  The  error  state  vector  holds  fifteen 
error  states  as  each  of  the  five  error  state  categories  are  propagated  in  three  directions 
(NED)  [28]  [30]. 

Sx(t)  =  F5x(t)  +  Bu(t)  +  Gw(t)  (3-17) 

Sx(t)  =  F  8x(t)  +  Gw  (t)  (3.18) 

5p  (t) 

5v(t) 

Sx(t)  =  ^(t)  (3.19) 

5a(t) 

5b(t) 

The  state  space  form  of  the  INS  error  state  motion  model  is  shown  in  Equation  3.20 
along  with  the  random  noise  inputs  [28]  [30].  The  position  states  do  not  contain  any 
input  noise  entering  into  them  due  to  fact  that  the  position  error  states  are  integrated 
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directly  from  the  velocity  error  states  which  do  contain  input  noise.  The  5v  and  ip 
error  states  both  are  random  walk  processes  [13]  [30].  The  5a  and  5b  error  states 
follow  a  FOGM  process  [13]  [30].  It  should  be  noted  that  the  system  noise  mapping 
matrix  G  is  located  along  the  second  line  of  Equation  3.20. 


5p(f) 

0 
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0 

5v(t) 
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-2Cnfi?  Ce 
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n>  J 

5p(f) 
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(3.20) 


The  input  noise  covariance  matrix  Q  (t)  for  the  error  state  motion  model  in  Equa¬ 
tion  3.21  has  input  noise  entering  into  the  last  four  error  state  categories  5v.  ip,  5a, 
and  5b  [30].  The  pre-dehned  error  noise  strength  values  for  these  states  used  in  INS 
error  state  simulation  are  presented  in  Table  3.1.  The  correlation  time  constant  r  for 
the  accelerometer  bias  5a  and  the  gyroscope  bias  5b  for  an  INS  are  3600  seconds  for 
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error  state  simulation  of  each  INS  grade  [13]. 


QW 


0  0  0  0  0 


0  0  0  0 

r 

o  0  0  0 

T 

0  0  0  0 

r 


0  0  0 


0 


^cr(5b(t) 

r 


(3.21) 


Table  3.1.  Noise  Strength  values  of  INS  Grades  for  use  in  error  simulation  [13] 


INS  Grade 

aSb(t)  C-f) 

aSa(t)  (32) 

^«)  ( f) 

/  m  \ 

ab(0  (s| ) 

Commercial  (Cloudcap  Crista) 

8.7  x  10“3 

1.96  x  10_1 

6.5  x  10"4 

4.3  x  10“3 

Tactical  (HG1700) 

4.8  x  10-6 

9.8  x  10"3 

8.7  x  10"5 

9.5  x  10-3 

Navigation  (HG  9900  -  H764G) 

7.2  x  10-9 

2.45  x  HP4 

5.8  x  10"7 

2.3  x  10-4 

Equation  3.23  is  the  b-frame  to  n-frame  DCM  used  in  the  error  state  motion  model 
[30].  The  inputs  to  the  DCM  are  vehicle  body  Euler  angle  rotations  roll  0,  pitch  9, 
and  yaw  0  which  are  measured  in  radians  as  shown  in  Equation  3.22.  The  DCM  can 
be  used  to  rotate  a  vector  from  the  b-frame  to  the  n-frame  which  is  in  meters  [28]  [30] . 
A  transpose  of  the  DCM  Cj]  will  result  in  the  reverse  DCM  C[], 


<P 


u>h 


0 


0 


(3.22) 


C 


n  _ 

b  — 


cos^cos^) 
sin('0)cos(0) 
— sin($) 


— sin('0)cos((/))  +  cos(^)sin(0)sin(</)) 
cos('0)cos((/))  +  sin(^)sin(0)sin(</)) 
cos(#)sin  (</>) 


sin('0)sin((/))  +  cos('0)sin(0)cos((/)) 
— cos('0)sin((/))  +  sin('0)sin(0)cos((/)) 
cos(#)cos(</>) 

(3.23) 
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Equation  3.24  is  the  DCM  that  rotates  position  information  from  the  e-frame  to  the 
n-frame  [28]  [30] .  The  DCM  takes  in  latitude  L  and  longitude  A  coordinates  in  radians. 

— sin(L)cos(A)  — sin(L)sin(A)  cos(L) 

C"  =  — sin(A)  cos  (L)  0  (3.24) 

— cos(L)cos(A)  — cos(L)sin(A)  —  sin(L) 

The  Earth’s  rotation  must  be  taken  into  account  when  computing  errors  for  an 
aircraft  traversing  over  the  Earth  [30].  Equation  3.25  represents  the  Earth’s  rotation 
with  respect  to  the  i- frame  [30].  In  the  INS  error  state  motion  model,  this  is  rotated 
into  the  n-frame  and  affects  the  tilt  error  state.  Equation  3.26  is  the  skew  symmetric 
matrix  of  Equation  3.25.  Equation  3.26  accounts  for  the  Earth’s  rotation  in  the 
velocity  error  states  [30].  The  Earth  rotation  rate  Q  used  in  this  research  effort  is 
7.292115  x  10“5^  (WGS-84)  [28]. 

(3.25) 

0  —Q  0 

n°e=  n  o  o  (3.26) 

0  0  0 

Gravity  over  the  surface  of  the  Earth  varies  slightly  based  upon  a  region’s  elevation 
[30].  To  account  for  this  phenomenon  within  the  position  error  component  of  the 
velocity  error  state,  a  gradient  value  D  is  used  as  a  unique  regional  approximation 
of  gravity  [30].  D  in  Equation  3.27  represents  a  gravitational  constant  unique  to  a 
specific  region  based  around  the  aircraft’s  current  position  with  respect  to  the  e-frame 
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pe  and  the  skew  symmetric  of  the  earth’s  rotation  vector  from  Equation  3.26  [30]. 
The  gravitational  constant  g  is  9.8  ^  [31]. 


D  =  V9|pl 


3p'(pe)T 


I 


(«L) 


2 


(3.27) 


Accelerometers  within  INS  units  output  specific  force  measurements  which  in  this 
research  effort  are  within  vector  fb  with  respect  to  the  b-frame  as  shown  in  Equation 
3.28  [30].  For  the  velocity  error  states,  the  specific  force  measurements  are  taken  into 
account  through  rotation  into  the  n- frame  as  seen  in  Equation  3.29  [30].  For  this 
research  effort,  the  simulated  aircraft  trajectory  is  a  wings  level  flight. 


fi  \ 

fn  =  C£fb 


(3.28) 


(3.29) 


3.5  Error  State  Motion  Model  Validation 

To  verify  the  INS  error  state  motion  model,  it  is  incumbent  to  achieve  a  match 
between  the  simulated  error  ensemble  statistics  of  the  INS  error  model  and  the  INS 
error  theoretical  statistics  [20].  The  theoretical  or  analytical  statistics  are  the  INS 
grade  dependent  theoretical  results  that  would  be  expected  of  INS  errors  propagated 
over  time.  The  theoretical  statistics  create  a  template  that  shows  expected  INS  errant 
behavior  over  time.  The  form  of  the  theoretical  statistics  can  be  inferred  due  to  the 
Gaussian  nature  of  the  model’s  distribution  where  the  mean  is  zero  and  the  standard 
deviation  follows  the  a  values  listed  in  Table  3.1  [20]. 

The  INS  error  state  motion  model’s  simulated  error  statistics  are  collected  from 
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propagating  multiple  runs  of  the  error  state  motion  model  over  time  in  Equation 
3.20  and  then  determined  by  the  ensemble  mean  and  standard  deviation  [5].  If  the 
ensemble  statistics  match  the  form  of  the  theoretical  statistics,  the  model  can  be 
considered  sound  as  a  means  to  simulate  INS  errors  during  a  vehicle  trajectory  [20]. 
To  prove  the  validity  of  the  error  state  motion  model  in  Equation  3.20,  a  Monte 
Carlo  simulation  of  10,000  runs  was  conducted  for  all  three  grades  of  INS  based  on 
the  simulation  parameters  listed  in  Table  3.2.  Figures  3. 3-3. 5  show  the  ensemble 
and  theoretical  statistics  matching  for  the  position  state  in  three  directions  for  a 
propagation  of  ten  minutes.  It  is  clear,  the  errors  decrease  exponentially  with  the 
change  in  the  quality  of  INS  grades  from  Commercial  to  Navigation  Grade. 

Table  3.2.  Simulation  Parameter  for  Error  State  Motion  Model  Validation  [13]. 


Monte  Carlo  Simulation  Parameters 


Propagation  Time  Step  (A t)k 

Is 

Aircraft  Initial  Position 

(0,0,0)  m 

Aircraft  Initial  Velocity 

(200,0,0)  f 

Aircraft  Initial  Acceleration 

(0,0,0)  f2 
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D  Error  (m)  E  Error  (m)  N  Error  (m) 


x  10s 


Time  (min) 


x  104 


Time  (min) 


Figure  3.3.  Monte  Carlo  Simulation  of  10,000  runs  for  Commercial  Grade  INS  Position 
Error  for  the  n-frame  North,  East,  Down  axes  respectively. 
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Figure  3.4.  Monte  Carlo  Simulation  of  10,000  runs  for  Tactical  Grade  INS  Position 
Error  for  the  n-frame  North,  East,  Down  axes  respectively. 
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Figure  3.5.  Monte  Carlo  Simulation  of  10,000  runs  for  Navigation  Grade  INS  Position 
Error  for  the  n-frame  North,  East,  Down  axes  respectively. 


3.6  Measurement  Model 

For  this  research  effort  INS  drift  errors  are  corrected  within  an  EKF  through  the 
use  of  lunar  feature  angle  measurements  obtained  in  the  n-frame  and  a  barometer. 

Lunar  Feature  Angle  Measurement  Model. 

Equation  3.30  shows  the  non-linear  measurement  model  containing  measurements 
zopt ,ki  the  measurement  function  hopt(-),  and  measurement  noise  vopt:/,:  which  is  a 
white  Gaussian  noise  process  [21].  Ropt,fc  is  the  measurement  noise  covariance  matrix. 
The  position  states  p?  are  the  input  to  the  measurement  function.  The  measurement 
function  hopt(-)  consists  of  two  external  angle  measurements:  azimuth  r;  and  elevation 
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£  described  in  Equations  3.33  and  3.34  respectively.  The  angles  are  calculated  based 
upon  the  relationship  of  the  simulated  vehicle  position  p /,.  to  simulated  lunar  features 
Tn  in  the  n- frame.  hopt(-)  is  constructed  to  allow  for  angles  to  be  calculated  based 
upon  multiple  features  for  batch  updating  denoted  by  index  1  to  n  with  discrete  times 
at  index  k. 

^opt ,k  h0pt(pfc)  T  (3.30) 

Ropt,fc  =  E[vopt,fcVoptA.]2nx2ri  (3.31) 


^opt(Pfc) 


Vi,k 

Vn,k 

€i,k 

£n,k 


J  2nxl 


(3.32) 


rj  is  determined  by  the  slope  of  the  vehicle  position  (pn,Pe,Pd)  in  relation  to  the 
feature  position  (T$,  T£,T£)  within  the  inverse  tangent  function  along  the  horizontal 
E-N  plane  as  shown  in  Equation  3.33  and  Figure  3.6.  The  positive  E-axis  is  07r 
rads.  £  is  determined  based  upon  the  slope  of  the  vehicle  position  in  relation  to  the 
feature  position  with  the  rise  along  the  vertical  D-axis  and  the  run  calculated  as  a 
line  projection  along  the  horizontal  E-N  plane  as  shown  in  Equation  3.34  and  Figure 


3.7. 


Vk  =  tan  1 


f  ?N  ~  PN,A 

Ve“  Pe ,k) 


(3.33) 


=  tan  1 


1D 


Pd,/c 


yuE"  -  pE,t)2  +  ps 


(3.34) 
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Figure  3.6.  3D  Angle:  E-N  Plane. 
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Figure  3.7.  3D  Angle:  D-axis  and  E-N  Plane  Projection. 
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For  EKF  estimation,  simply  inputing  the  non-linear  measurement  model  into  the 
filter  will  not  work  due  to  the  linear  nature  of  the  Kalman  gain  and  covariance  update 
equations  [21].  The  angle  measurements  must  be  linearly  approximated  along  the 
three  orthogonal  axes  that  comprise  the  n-frame.  A  first-order  Taylor  series  expansion 
called  a  Jacobian  is  used  to  linearly  approximate  the  non-linear  measurement  function 
as  seen  in  Equation  3.36  [21].  Within  the  EKF  algorithm,  the  Jacobian  is  determined 
for  the  estimated  angles  fj  and  £  which  have  been  calculated  based  upon  estimated 
position  states  pk  shown  in  Equation  3.35. 


P  k  =  PiNS,fc  +  Spk  (3.35) 

Hopt,fc  =  w?-hopt(pA;)  (3.36) 

op  k 

Hopt;fc  represents  the  measurement  matrix  comprised  of  Jacobians  required  for  the 
EKF  update  equations,  rj  and  £  only  give  observability  along  the  position  states. 
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(3.37) 
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For  fj,  the  Jacobians  along  the  north  and  east  position  states  are  represented  below. 


H 


°pt,pN,fc 


d  ,  -(Tj  -  pE,fc) 

df>N,k^1,k  (Tg  —  pE,fc)2  +  (Tn  —  Pn,a;)2 


(3.38) 


H 


opt-PE.fe 


d  ,  _  _ (Tn  —  PN,fc) _ 

dpE,kVl,k  (Te  —  PE,fc)2  +  (?N  —  PN,fc)2 


(3.39) 
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For  £,  the  Jacobians  along  the  north,  east,  and  down  position  states  are  represented 
below. 


Hopt,pN,fc  - 


_  (PN,fc  M  PD,fc) _ 

(Tg  —  PE,fc)2  +  (^N  —  PN ,k)2  ((^E  _  PE,fc)2  +  (2~n  —  PN,fc)2  +  (^D  —  PD,fc)^ 

(3.40) 


lopt,pE  ( 


-£l  m  — 


_ (pe^  Tg)(rg  pp,fe) _ 

(Te  —  Pe,^)2  +  (^N  —  PN,fc)2  m  —  Pe ,k)2  +  (^N  —  PN,fc)2  +  (^D  —  Pd ,kY 

(3.41) 


H°pt,pDifc 

For  this  research  effort,  the  standard  deviation  <jopt  for  lunar  feature  angle  mea¬ 
surement  noise  vopt  /,.  is  set  as  the  angular  resolution  A6q  of  features  in  an  image 
constrained  by  lens  diffraction  as  discussed  in  Section  2.5.  To  design  the  camera  de¬ 
tector  to  the  diffraction  limit,  an  assumption  was  made  based  on  a  common  optical 
design  that  sampling  should  be  equal  to  the  diffraction  limited  spot  size.  With  this 
stipulation,  the  Airy  disk  diameter  2qi  can  be  assumed  to  equal  the  pixel  size  At  for 
the  diffraction  limited  system  used  in  this  research  effort  as  shown  in  Equation  3.43 
and  Figure  3.8.  This  assumption  places  the  angular  resolution  within  ±1  pixel.  To 
obtain  cropt,  Equation  3.44  presents  the  relationship  between  the  pixel  size  and  the 
focal  length  used  to  derive  A9q. 

2 qx  «  At  (3.43) 


'6 ,k  ~ 


\j  (^E  PE,fc)2  +  (2~n  Pn,/c)2 

(Tg  -  PE,fc)2  +  (T£  -  Pn ,k)2  +  (T£  -  p D,fe) 


(3.42) 


(3.44) 
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Pixel  Squares 


M 


Figure  3.8.  Diffraction  Limited  Spot  diameter  assumed  equal  to  Pixel  Size. 

Barometer  Measurement  Model. 

The  down  position  state  is  completely  observed  by  a  barometer  which  measures 
the  aircraft’s  altitude  as  shown  in  Equations  3.45  and  3.46  [5].  For  this  research  effort, 
the  barometer  measurements  are  brought  into  the  EKF  every  60  seconds. 


^baro ,k  Hbaro,/c-X-fc  T  Vbaro,/c 


(3.45) 


H 


baro,fc 


0  0  10...  0 


(3.46) 


L  J  1x15 

The  barometer  measurement  noise  Vbaro,fc  is  a  white  Gaussian  process.  The  standard 
deviation  (j|)aro  for  the  barometer  measurement  noise  is  set  at  ±5m.  The  measurement 
covariance  matrix  Rbaro,fc  is  represented  below. 


Rbaro,fc  —  E  [vba.ro. A:  vba.ro./c] 


15x15 


(3.47) 
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3.7  Feature  Design  and  Implementation 


To  simulate  features  derived  from  the  lunar  surface,  randomized  azimuth  and  ele¬ 
vation  angles  are  first  generated  in  vectors  0feat  and  t/>feat  using  the  Equations  3.48  and 
3.49  respectively.  These  equations  represent  a  common  practice  in  placing  points  at 
random  locations  along  a  spherical  surface  based  on  the  spherical  coordinate  system 
equations.  mfeat  and  nfeat  are  vectors  containing  random  numbers  with  range  between 
0  and  1.  The  angles  are  then  processed  through  the  MATLAB®  function  sph2cart 
with  the  Moon’s  radius  to  gain  cartesian  coordinate  feature  positions  (T£,  Ty,  X|) 
in  the  ECEF  frame.  The  cartesian  coordinates  are  then  translated  to  the  Moon’s 
position  in  orbit  based  upon  the  orbital  parameters  in  Section  2.6  providing  an  ap¬ 
proximation  of  lunar  surface  features.  Figures  3.9-3.11  present  1000  simulated  lunar 
features  as  green  x’s  and  the  Moon  as  a  red  spheroid. 


@  feat  Hlfeat7r 


(3.48) 


^feat  =  sin  *(-l  +  2nfeat) 


(3.49) 


^  1.855  . 


g  1 .835  . 
0-  1.83. 


^  1 .825  . 


X  Position  (km) 


Figure  3.9.  ECEF  Frame:  Moon  and  1000  Random  Lunar  Features. 
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Y  Position  (km) 

Figure  3.10.  ECEF  Frame:  Y-Z  Plane.  Moon  and  1000  Random  Lunar  Features. 
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Figure  3.11.  ECEF  Frame:  X-Z  Plane.  Moon  and  1000  Random  Lunar  Features. 

3.8  EKF  Implementation 

This  section  provides  further  details  about  the  system  algorithm  setup  and  EKF 
implementation.  To  begin,  orientation  of  the  aircraft’s  truth  trajectory  and  the 
INS  trajectory  xjns ,k  is  needed  in  the  context  of  traveling  along  the  Earth’s  surface. 
For  this  reason,  the  algorithm  requires  a  starting  position  in  geodetic  coordinates: 
latitude,  longitude,  and  altitude  (L,  A,/i).  The  initial  geodetic  coordinate  incorpora¬ 
tion  into  the  DCM  helps  to  orient  the  fixed  n-frame  of  the  aircraft’s  trajectory  to 
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the  Earth’s  surface.  Outside  of  EKF  error  state  estimation  propagation,  the  initial 
geodetic  coordinate  is  solely  used  in  the  C®  DCM  to  orient  the  n-frame  due  to  the 
small  aircraft  trajectory  modeled.  The  resulting  posture,  places  the  fixed  n-frame’s  N 
and  E  axes  as  a  plane  tangent  to  the  Earth’s  surface  resting  atop  the  starting  geodetic 
coordinate  where  the  D-axis  points  toward  the  Earth’s  center. 

An  overall  system  representation  for  Lunar  Vision  Aided  Navigation  with  Inertial 
Navigation  System  is  shown  in  Figure  3.12.  The  figure  shows  the  INS  trajectory 
xiNS,fc  is  rotated  from  the  n-frame  into  geodetic  coordinates  and  input  into  the  EKF 
to  estimate  INS  error  states  hx*,  and  a  trajectory  estimate  x*  in  the  n-frame.  The 
simulated  lunar  features  are  in  the  e- frame  to  best  represent  their  orientation  to  Earth. 
The  lunar  features  are  randomized  along  the  Moon’s  near  side  surface  facing  Earth. 
To  perform  the  r)  and  £  angle  calculations  and  the  Jacobians  in  Equations  3.38-3.42 
within  the  EKF,  the  features  and  the  trajectory  estimate  need  to  both  be  within 
the  same  frame.  To  accomplish  this,  the  lunar  features  are  rotated  into  the  n-frame 
based  upon  the  initial  geodetic  coordinate.  Once  the  EKF  estimation  is  complete, 
the  estimated  trajectory  can  be  rotated  into  the  e- frame  to  ECEF  coordinates  or 
converted  to  geodetic  coordinates  to  gain  knowledge  of  the  vehicle  trajectory  along 
the  surface  of  the  Earth. 
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♦Lunar  Feature  Angles 
♦Barometer 


Figure  3.12.  Lunar  Vision  Aided  Navigation  with  Inerital  Navigation  System  [5]. 


It  should  be  noted  that  based  on  the  nature  of  the  Jacobians  described  in  Section 
3.6,  an  angle  obtained  by  a  vehicle  observing  a  feature  at  infinity  within  the  same 
reference  frame  will  appear  to  be  constant  at  90°.  The  lack  of  change  in  the  observed 
angle  will  cause  the  Jacobians  to  converge  to  zero  resulting  in  an  observability  loss  of 
the  vehicle’s  position  error  states.  In  this  research  effort,  however,  the  fixed  n-frame 
along  the  Earth’s  surface  and  the  lunar  feature  finite  distances  allows  for  observability 
in  the  position  states.  The  down  position  state  is  observed  fully  by  a  barometer. 
Thus,  INS  drift  error  is  corrected  in  all  position  states.  The  loss  of  position  state 
observability  by  a  decrease  in  Jacobians  based  on  feature  distance  will  be  presented 
Chapter  4.  Also  in  Chapter  4,  a  singularity  in  the  spherical  coordinate  system  where 
this  concept  does  not  hold  true  will  be  shown  as  well. 


3.9  Summary 

Chapter  3  provided  the  research  methodology  to  obtain  the  results  in  Chapter  4. 
The  truth  model,  error  state  motion  model,  measurement  model,  and  lunar  feature 
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model  were  presented.  The  error  state  motion  model  validation  was  shown,  and  the 
EKF  implementation  was  discussed. 
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IV.  Results  and  Analysis  of  the  Theoretical  Limits  of  Lunar 
Vision  Aided  Navigation  with  INS 

4.1  Introduction 

This  chapter  presents  the  results  from  the  research  methodology  detailed  in  Chap¬ 
ter  3.  The  purpose  of  this  research  was  to  simulate  and  analyze  INS  drift  correction 
by  lunar  feature  angle  measurements.  Various  scenarios  were  considered  in  order  to 
better  understand  the  impact  that  lunar  feature  aided  navigation  would  theoretically 
have  on  INS  drift  correction.  The  affect  of  lunar  feature  aided  navigation  on  different 
INS  grades  is  shown  in  Section  4.2.  Measurement  noise  effects  from  camera  specifi¬ 
cations  on  INS  drift  correction  by  lunar  feature  angle  measurements  are  presented  in 
Section  4.3.  Feature  quantity  effects  are  shown  in  Section  4.4.  Section  4.5  covers  the 
effect  of  lunar  features  at  increasing  distances  and  different  lunar  orbit  geometries. 
Section  4.6  presents  the  chapter  summary. 

It  must  be  stated  that  for  the  scenarios  in  Chapter  4,  lunar  vision  aided  navigation 
with  INS  was  simulated  under  an  idealized  set  of  parameters  where  measurement  noise 
effects  are  only  considered.  Therefore,  the  results  in  Chapter  4  are  not  indicative  of 
actual  performance  by  lunar  vision  aided  navigation  with  INS  but  present  analysis  on 
the  trade  spaces  resulting  from  changes  to  INS  grade,  camera  image  array  pixel  sizes, 
camera  lens  focal  length,  lunar  feature  quantity,  and  lunar  feature  distance  and  ge¬ 
ometry.  The  many  assumptions  which  allow  the  aforementioned  model  simplification 
are  detailed  in  Section  1.4. 

For  the  following  scenarios,  unless  stated  otherwise,  the  parameters  in  Table  4.1 
are  set  for  the  simulated  trajectory  runs.  The  default  camera  measurement  noise 
standard  deviation  <7opt  is  set  at  O.Ollmrads  based  on  the  angular  resolution  Adq  of  a 
Prosilica  GE1660C  CCD  camera  with  a  253mm  focal  length  lens.  A  barometer  aids 
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the  down  position  state  with  a  measurement  noise  standard  deviation  a^aro  set  at 
±5m. 


Table  4.1.  Default  Case  Trajectory  Parameters 


Samples 

Time 

Start  Pos*  ( LXh ) 

INS 

Prop*  Time 

■#-  of  Features 

Camera  Update 

Baro  Update 

30 

lOmin 

(0°,0°,0m) 

Tactical 

Is 

100 

4s 

60s 

*  Prop  -  Propagation,  Pos  -  Position 


Figures  4. 1-4.4  present  the  default  orientation  of  the  Moon  to  the  Earth  in  the 
ECEF  frame.  The  blue  sphere  represents  the  Earth  and  the  dark  red  smaller  sphere 
represents  the  Moon.  In  the  default  orientation,  the  Moon  is  at  its  highest  inclination 
to  the  Earth’s  equator  at  28.58°.  The  Moon  lies  at  an  offset  above  the  X-axis.  In 
reference  to  the  true  trajectory  in  the  n-frame  along  the  Earth’s  surface,  the  course 
of  the  aircraft  travels  in  a  circle  with  radius  20km  along  the  N-E  plane  for  every 
sample  in  the  default  case.  The  height  of  the  aircraft  varies  based  upon  the  FOGM 
acceleration  model  detailed  in  Section  3.3. 


Figure  4.1.  ECEF  Frame.  Earth  and  Moon  with  fixed  local  level  n-frame  (Not  to 
Scale). 
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n-frame 


Figure  4.2.  ECEF  Frame:  X-Z  Plane.  Earth  and  Moon  with  fixed  local  level  n-frame 
(Not  to  Scale). 
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Figure  4.3.  ECEF  Frame:  X-Y  Plane.  Earth  and  Moon  with  fixed  local  level  n-frame 
(Not  to  Scale). 
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Figure  4.4.  ECEF  Frame:  Y-Z  Plane.  Earth  and  Moon  with  fixed  local  level  n-frame 
(Not  to  Scale). 
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4.2  INS  Grade  Testing 


This  section  presents  INS  drift  correction  aided  by  lunar  feature  angles  against 
varying  grades  of  INS.  All  realizations  were  ran  using  the  parameters  in  Table  4.1 
with  100  features,  and  the  camera  measurement  noise  standard  deviation  was  kept 
at  O.Ollmrads.  Figures  4. 5-4. 7  show  significant  decrease  between  the  position  state 
error  for  a  stand  alone  INS  when  compared  to  the  position  state  errors  from  an  INS 
aided  by  lunar  vision  angle  measurements.  In  addition,  the  position  state  errors  show 
significant  decrease  with  each  change  in  INS  grade  when  comparing  Figures  4.5b,  4.6b, 
and  4.7b.  For  the  Commercial  Grade  INS  in  Figure  4.5b,  the  north  position  state 
error  variance  is  at  ±58. 96m  and  the  east  position  state  error  variance  is  at  ±25. 85m. 
The  down  position  state  error  variance  is  sharply  high  at  ±77. 62m  even  with  the 
barometer  having  full  observability  of  the  down  position  state.  For  the  Navigation 
Grade  INS  in  Figure  4.7b,  the  error  significantly  decreases  for  the  position  states 
with  error  variances  at  ±2. 946m  for  the  north  position  state,  ±0. 9805m  for  the  east 
position  state,  and  ±2. 328m  for  down  position  state. 


Time  (min) 


(a)  No  Lunar  Vision  Aiding  (b)  Lunar  Vision  Aiding 

Figure  4.5.  Commercial  Grade  INS  with  100  Features.  Position  State  Error  for  the 
n-frame  North,  East,  Down  States  respectively. 
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Figure  4.6.  Tactical  Grade  INS  with  100  Features.  Position  State  Error  for  the  n- frame 
North,  East,  Down  States  respectively. 


Time  (min)  Time  (min) 


(a)  No  Lunar  Vision  Aiding 


(b)  Lunar  Vision  Aiding 


Figure  4.7.  Navigation  Grade  INS  with  100  Features.  Position  State  Error  for  the 
n-frame  North,  East,  Down  States  respectively. 


4.3  Camera  Specifications  Testing 

The  effect  of  camera  measurement  noise  vopt^  on  an  EKF  estimating  INS  errors 
aided  by  lunar  feature  angle  measurements  can  be  predicted.  A  relationship  exists 
within  an  EKF  whereby  the  measurement  noise  amount  can  cause  the  EKF  to  rely 
more  on  the  current  error  state  estimate  or  rather  on  the  correction  to  the  current 


error  state  estimate  brought  about  by  the  incorporated  measurement  z/,.  as  shown  in 
Section  2.4. 

An  EKF  incorporating  measurements  with  high  noise  will  cause  a  low  Kalman 
gain.  The  low  Kalman  gain  will  decrease  the  amount  of  correction  by  the  measure¬ 
ments  causing  the  EKF  to  rely  more  on  the  current  error  state  estimate.  For  this 
research  effort,  Figures  3. 3-3. 5  in  Section  3.5  present  INS  error  states  without  correc¬ 
tion.  For  lower  measurement  noises,  the  error  state  estimates  will  be  much  closer  to 
the  true  error  states  relating  to  higher  accuracy  in  estimation.  Lower  measurement 
noises  allow  for  a  higher  quality  estimate  correction  and  a  Kalman  gain  which  allows 
for  the  full  proportion  of  the  near  accurate  correction  to  update  the  current  state 
estimate. 

For  this  research  effort,  the  standard  deviation  aopt  of  the  camera  measurement 
noise  that  is  added  to  the  observed  angle  measurements  is  obtained  from  the  angular 
resolution  A 0q  of  a  diffraction  limited  system  which  is  a  fixed  camera  mounted  to  an 
aircraft  body  frame.  It  was  assumed  that  the  diameter  of  a  diffraction  limited  spot 
2qi  is  set  equal  to  the  pixel  size  A£  of  the  camera’s  image  array  so  as  to  constrain  the 
image  uncertainty  to  ±1  pixels.  Section  3.6  presented  a  relationship  connecting  pixel 
size  At  to  the  angular  resolution  Adq.  Table  4.2  presents  a  sweep  of  pixel  sizes  and 
the  coinciding  angular  resolution  to  test  INS  drift  correction  by  lunar  feature-aided 
navigation  against.  The  focal  length  is  set  at  253mm.  To  place  the  pixel  sizes  into 


Table  4.2.  Measurement  Noise 


Measurement  Noise  Level 

A£  ( /jm) 

v0pt,/c  (mrads) 

1 

1000 

2 

2 

100 

0.2 

3 

10 

0.02 

4 

1 

0.002 

perspective,  camera  image  sensor  pixel  sizes  were  obtained  from  commercial  camera 
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specifications.  For  the  SCI  High  Accuracy  Star  Tracker  2  image  sensor,  the  pixel  size 
was  18/jin.  This  is  near  the  pixel  size  for  the  Prosilica  GE1660C  camera  which  was 
5.5/tm. 

From  the  previous  assertions,  Figures  4.8-4.11  present  the  position  state  error  in 
the  n-frame  with  corrections  from  the  angle  measurements.  An  accompanying  single 
run  of  an  estimated  trajectory  tracking  a  true  trajectory  is  presented  as  well  in  the 
right  subfigure  within  Figures  4.8-4.11.  All  of  the  trajectories  shown  are  based  upon 
the  parameters  listed  in  Table  4.1. 

With  each  decrease  in  pixel  size  and  measurement  noise,  the  position  state  error 
decreases  for  the  north  and  east  position  states.  In  Figure  4.8a,  north  position  state 
error  variance  is  ±396m  and  east  position  state  error  variance  is  ±2 16. 5m  based 
on  a  measurement  noise  of  2mrads  and  pixel  size  at  1000/im.  In  contrast,  Figure 
4.11a  presents  north  position  error  variance  at  ±8. 33m  and  ±1.54m  for  the  east 
position  state  for  a  measurement  noise  of  0.002mrads  and  pixel  size  l//rn.  Based  upon 
this  information,  the  results  confirm  the  predicted  increase  in  state  error  estimation 
accuracy  based  upon  decreasing  pixel  size  within  this  system  model. 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.8.  Pixel  Size:  1000pm  and  Measurement  Noise:  2mrads. 
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(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.9.  Pixel  Size:  100/im  and  Measurement  Noise:  0.2mrads. 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.10.  Pixel  Size:  10/im  and  Measurement  Noise:  0.02mrads. 
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(b)  Estimated  vs.  True  Trajectory 


(a)  Position  State  Error  for  the  n-frame  North, 

East,  Down  axes  respectively 

Figure  4.11.  Pixel  Size:  1pm  and  Measurement  Noise:  0.002mrads. 


Commercial  Camera  Lenses. 

For  the  below  scenarios,  the  parameters  from  a  commercial  camera  and  various 
lenses  were  obtained  to  observe  their  behavior  in  INS  drift  correction  by  lunar  vision 
aided  navigation.  The  trajectory  runs  followed  the  parameters  listed  in  Table  4.1. 
The  camera  used  in  the  simulations  was  a  Prosilica  GE1660C  CCD  Camera  with  a 
pixel  size  listed  in  the  maintenance  guide  as  5.5//mx5.5//m.  From  the  relationship 
established  in  Section  3.6  between  lens  focal  length  /  and  the  angular  measurement 
noise  standard  deviation  <ropt,  a  pattern  can  be  shown  that  as  /  increases  the  cropt 
decreases.  This  results  in  a  better  trajectory  estimation  within  this  model.  Table  4.3 
presents  the  lens  parameters  for  lenses  that  are  compatible  with  the  c- mount  2/3” 
optical  format  of  the  Prosilica  GE1660C  CCD  camera. 

The  image  array  is  rectangular  and  has  an  associated  rectangular  field  of  view 
(FOV)  with  width  angle  aw  and  height  angle  a^.  Based  on  the  field  of  view  equations 
in  Section  2.5,  the  planar  area  of  the  observed  world  scene  by  the  camera  can  be 
calculated  as  shown  in  the  last  column  listing  the  observed  world  scene  width  wscene 
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and  height  hscene  parameters.  The  first  two  lens  FOV  areas  listed  can  completely 
view  the  Moon’s  planar  projection  within  their  width  and  height  parameters.  The 
Moon’s  planar  projection  in  the  world  scene  viewed  by  the  camera  appears  as  a 
circle  with  diameter  3475km.  The  FOV  area  of  the  Cinetel  III  Telephoto  Lens  falls 
short  of  completely  viewing  the  planar  projection  of  the  Moon.  However,  the  zoom 
capabilities  of  the  Cinetel  III  Telephoto  Lens  with  focal  length  1000mm  allows  for  a 
much  closer  view  of  the  Moon’s  features  leading  to  a  smaller  measurement  noise  vopt,A- 
and  a  tighter  estimate. 


Table  4.3.  Commercial  Lens  Parameters 


Camera  Lens 

Focal  Length 
(mm) 

^opt 

f-number 

Lens  Diameter 
(mm) 

«h(°) 

«w(°) 

FOV  Area 
(km) 

Tamron  23FM50SP 

50 

0.055 

f/1.4 

35.7 

7.6 

10.06 

67,654.4x50,740.8 

Fujinon  22x  Telephoto  Zoom  Lens 

253 

0.011 

f/1.6 

158.12 

1.5 

2 

13,370.4x10,027.83 

Cinetel  III  Telephoto  Lens 

1000 

0.003 

f/5.6 

178.57 

0.4 

0.5 

3,382.72x2,537.04 

Figures  4.12-4.14  present  the  position  state  error  based  upon  the  Prosilica  GE1660C 
CCD  Camera  with  the  lenses  listed  in  Table  4.3  and  the  associated  measurement  noise 
a  listed  in  the  respective  figure  captions.  In  Figure  4.12a  the  north  position  state  er¬ 
ror  variance  is  at  ±28. 67m  and  the  east  position  state  error  variance  is  at  ±11. 25m. 
With  the  decrease  in  measurement  noise  and  the  increase  in  focal  length  to  obtain  a 
closer  view  of  the  Moon’s  features,  Figure  4.14a  shows  the  north  position  state  error 
variance  decreasing  to  ±6. 42m  and  the  east  position  state  error  variance  decreasing 
to  ±2. 347m. 
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(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.12.  Focal  Length:  50mm  and  Measurement  Noise  a:  0.055mrads. 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.13.  Focal  Length:  253mm  and  Measurement  Noise  a:  O.Ollmrads. 
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(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.14.  Focal  Length:  1000mm  and  Measurement  Noise  a:  0.003mrads. 


4.4  Feature  Quantity  Testing 

The  following  section  explores  the  accuracy  of  the  trajectory  estimation  predicated 
on  the  amount  of  lunar  features  observed.  Based  upon  the  EKF  algorithm,  as  a 
batch  angle  measurement  update  increases  in  the  number  of  angles  measured,  the 
measurement  model  z*,  increases  and  allows  for  multiple  angles  to  be  used  in  creating 
the  estimate  correction  and  provides  multiple  position  state  observabilities  in  the 
measurement  matrix  H&.  The  combination  of  the  Kalman  gain  and  the  residual  built 
upon  multiple  measured  angles  provides  an  increasingly  accurate  correction  to  the 
estimated  error  states  resulting  in  a  more  accurate  track  of  the  true  trajectory  by  the 
estimated  trajectory. 

From  Figures  4.15-4.22,  a  pattern  can  be  seen  in  the  increase  in  accuracy  as  the 
amount  of  lunar  features  increases.  For  the  single  feature  in  Figure  4.15,  it  can  be 
viewed  as  the  green  x  along  the  lower  portion  of  the  red  sphere  facing  the  ECEF  Y- 
axis.  In  Figure  4.16a,  the  north  position  state  error  variance  at  1  feature  is  ±273. 9m 
and  the  east  position  state  error  variance  is  ±114m.  The  error  variance  significantly 
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decreases  to  ±7. 435m  for  the  north  position  state  error  variance,  and  the  east  position 
state  error  variance  decreases  to  ±1.326m  as  shown  in  Figure  4.22a.  The  decrease 
in  the  error  variance  for  the  position  states  confirms  the  previous  prediction  that  the 
increase  in  features  improves  the  EKF  error  estimate.  The  trajectory  track  in  Figure 
4.22b  shows  what  appears  to  be  a  loose  track  of  the  true  trajectory  when  compared  to 
the  trajectory  estimation  with  less  features  in  Figure  4.18b.  However,  the  associated 
position  state  error  plots  in  Figure  4.22a  for  1000  features  confirms  the  substantial 
decrease  in  position  state  error  with  the  increase  to  1000  features.  The  drift  in  Figure 
4.22b  can  then  be  attributed  to  only  the  down  position  state  which  is  observed  by 
a  barometer  which  provides  the  same  down  position  error  correction  for  all  of  the 
feature  quantity  testing  scenarios. 


5 

x  10 


343  X  Position  (km) 

Figure  4.15.  ECEF  Frame.  Moon  and  1  Random  Feature. 
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Down  Error  (m)  East  Error  (m)  North  Error  (m) 


(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.16.  Feature  Quantity:  1. 
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X  10 


343  X  Position  (km) 

Figure  4.17.  ECEF  Frame:  Moon  and  10  Random  Features. 
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Down  Error  (m)  East  Error  (m)  North  Error  (m) 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.18.  Feature  Quantity:  10. 


x  10* 


Figure  4.19.  ECEF  Frame.  Moon  and  100  Random  Features. 
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Down  Error  (m)  East  Error  (m)  North  Error  (m) 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.20.  Feature  Quantity:  100. 


Figure  4.21.  ECEF  Frame.  Moon  and  1000  Random  Features. 
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(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.22.  Feature  Quantity:  1000. 


4.5  Feature  Distance  and  Geometry  Testing 

The  following  set  of  scenarios  that  were  conducted  are  in  reference  to  performing 
INS  drift  correction  aided  by  lunar  features  placed  at  various  distances  and  lunar 
geometries  to  show  the  effect  upon  the  n-frame  position  state  estimates.  All  trajectory 
runs  are  based  on  parameters  established  in  Table  4.1. 

Feature  Distance  Testing. 

In  this  research  effort,  an  EKF  incorporates  lunar  feature  angle  measurements  as 
a  correction  to  INS  errors.  As  discussed  in  Section  3.8,  feature  distances  can  effect 
the  amount  of  information  the  measurements  provide  about  the  states  to  the  EKF 
for  error  state  correction.  This  is  due  to  the  fact  that  Jacobians  based  on  change  in 
angle  with  respect  to  traveled  distances  are  used  to  provide  linear  observability  of  the 
position  states.  Measured  angles  of  a  feature  at  infinity  within  the  same  reference 
frame  all  converge  to  90°.  The  differences  between  the  measured  angles  of  a  feature 
at  infinity  become  indiscernible.  The  resulting  calculated  Jacobians  converge  to  zero 
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providing  no  observability  of  the  position  states  in  the  measurement  matrix  Hopt  /,:  and 
no  error  state  correction.  Feature  distances  and  orientation  causing  unobservability 
of  each  position  state  distinctly  can  occur  as  well  when  features  are  placed  straight 
along  an  n-frame  axis  at  a  great  distance.  For  features  located  at  a  far  distance  along 
an  n-frame  axis,  no  change  in  angle  occurs  with  respect  to  the  distance  traveled  along 
the  said  axis  resulting  in  a  lack  of  observability  and  correction  along  the  particular 
position  state. 

The  following  scenarios  seek  to  show  that  though  features  at  infinity  cause  unob- 
servability,  the  features  located  at  the  Moon’s  distance  still  provide  usable  Jacobian 
values.  Also,  the  effect  of  feature  distances  located  along  an  n-frame  axis  causing 
unobservability  for  the  associated  position  state  distinctly  will  be  shown.  For  the 
following  simulations,  the  Sun’s  distance  is  considered  infinity  for  simulation.  The 
ECEF  frame  axes  are  denoted  by  X  ,Y,  Z  and  the  fixed  local  level  n-frame  axes  are 
N,  E,  D  as  discussed  in  Section  2.3. 

n-frame  N-axis  Distance  Testing. 

Figures  4.24-4.29  show  the  results  of  INS  drift  correction  aided  by  lunar  features 
placed  along  the  north  axis  as  shown  in  Figure  4.23.  The  orientation  of  the  features 
in  Figure  4.23  in  no  way  represents  the  orientation  of  the  Moon’s  orbit  around  the 
Earth.  The  orientation  of  the  features  along  the  ECEF  Z-axis  is  set  up  to  provide 
insight  into  INS  drift  correction  behavior  aided  by  lunar  features  located  along  the 
N-axis  of  the  n-frame.  Figure  4.24a  confirms  the  assertion  that  the  features  directly 
along  a  single  axis  will  cause  a  decrease  in  observability  of  the  position  state  along 
that  axis.  The  position  state  error  variance  in  Figure  4.24a  shows  a  much  higher  error 
variance  for  the  north  position  state  than  the  east  position  state.  This  is  due  to  the 
lack  of  angular  changes  of  the  aircraft’s  north  position  with  respect  to  the  features 
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along  the  N-axis  which  causes  a  low  observability  of  the  north  position  state  by  the 
Jacobian.  The  pattern  is  repeated  in  Figure  4.25a  for  the  features  being  moved  along 
the  positive  ECEF  Z-axis  to  the  Sun’s  distance.  Also,  at  this  distance,  the  north 
position  state  error  variance  grows  to  the  form  of  the  uncorrected  drift  error  variance 
for  the  north  position  state  error  simulated  in  Section  3.5  for  a  Tactical  Grade  INS. 
The  east  position  state  error  also  begins  to  diverge  back  to  the  original  uncorrected 
drift  shown  in  Section  3.5  for  a  Tactical  Grade  INS  as  well. 
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Figure  4.23.  ECEF  Frame.  Features  along  n-frame  N-axis  and  ECEF  Z-axis  (Not  to 
Scale). 
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(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.24.  Features  Moon  Distance  along  N-axis  and  ECEF  Z-axis. 


(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 

Figure  4.25.  Features  at  Sun  Distance  along  N-axis  and  ECEF  Z-axis. 


The  ensemble  statistics  for  the  azimuth  angle  Jacobians  in  Figures  4.26-4.27  also 
show  the  decrease  in  observability  by  order  of  magnitude  for  the  north  position  error 
state  with  the  increase  in  distance  along  the  N-axis.  For  the  elevation  angle  Jacobian 
ensemble  statistics  in  Figures  4.28-4.29,  all  three  position  error  states  loose  observ¬ 
ability  by  an  order  of  magnitude  with  the  change  in  distance  to  the  Sun’s  distance. 
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(a)  Azimuth  Angle  Jacobians  for  the  n-frame  (b)  Scaled  View  of  Figure  4.26a 

North,  East,  Down  axes  respectively 

Figure  4.26.  Features  at  Moon  Distance  along  N-axis  and  ECEF  Z-axis. 
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(a)  Azimuth  Angle  Jacobians  for  the  n-frame  (b)  Scaled  View  of  Figure  4.27a 

North,  East,  Down  axes  respectively 

Figure  4.27.  Features  at  Sun  Distance  along  N-axis  and  ECEF  Z-axis. 
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(a)  Elevation  Angle  Jacobians  for  the  n- frame 


(b)  Scaled  View  of  Figure  4.28a 


North,  East,  Down  axes  respectively 


Figure  4.28.  Features  at  Moon  Distance  along  N-axis  and  ECEF  Z-axis. 


(a)  Elevation  Angle  Jacobians  for  the  n- frame 


(b)  Scaled  View  of  Figure  4.29a 


North,  East,  Down  axes  respectively 


Figure  4.29.  Features  at  Sun  Distance  along  N-axis  and  ECEF  Z-axis. 


n-frame  E-axis  Distance  Testing. 

For  the  next  two  scenarios  features  were  placed  along  the  n-frame  E-axis  at  the 
Moon’s  distance  for  the  first  scenario  and  then  at  the  Sun’s  distance  for  the  second 
scenario  to  observe  the  INS  drift  correction  with  lunar  features  directly  along  the  east 
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axis.  Figure  4.30  shows  the  orientation  of  the  features  with  respect  to  the  n-frame 
E-axis  which  coincides  with  the  ECEF  Y-axis  for  this  set  of  scenarios. 


Figure  4.30.  ECEF  Frame.  Features  along  n-frame  E-axis  and  ECEF  Y-axis  (Not  to 
Scale). 

Figures  4.31a  and  4.32a  confirm  the  claim  that  the  correction  to  INS  errors  by 
lunar  feature  angular  measurements  will  decrease  at  great  distances.  Both  position 
state  error  variances  in  Figure  4.32a  are  near  the  simulated  uncorrected  Tactical  Grade 
INS  position  state  error  variance  presented  in  Section  3.5.  Figure  4.31a  also  confirms 
the  relationship  of  the  lack  of  position  state  observability  with  features  at  the  lunar 
distance  along  the  E-axis  causing  a  larger  variance  for  the  the  east  position  state  error 
with  an  approximately  240m  greater  error  variance  compared  to  the  north  position 
state  error  variance.  The  azimuth  angle  Jacobian  ensemble  statistics  in  Figures  4.33 
and  4.34  confirm  a  decrease  in  observability  at  the  greater  distance  of  the  features 
along  the  E-axis  with  the  drop  in  observability  by  an  order  of  magnitude  for  the  east 
position  state  Jacobian.  The  elevation  Jacobians  in  Figures  4.35-4.36  also  show  an 
decrease  by  an  order  of  magnitude  for  the  three  position  state  errors. 
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(a)  Position  State  Errors  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 

Figure  4.31.  Features  at  Moon  Distance  along  E-axis  and  ECEF  Y-axis. 


(a)  Position  State  Errors  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.32.  Features  at  Sun  Distance  along  E-axis  and  ECEF  Y-axis. 
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(b)  Scaled  View  of  Figure  4.33a 


(a)  Azimuth  Angle  Jacobians  for  the  n- frame 
North,  East,  Down  axes  respectively 

Figure  4.33.  Features  at  Moon  Distance  along  E-axis  and  ECEF  Y-axis 
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(a)  Azimuth  Angle  Jacobians  for  the  n- frame 
North,  East,  Down  axes  respectively 
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(b)  Scaled  View  of  Figure  4.34a 


Figure  4.34.  Features  at  Sun  Distance  along  E-axis  and  ECEF  Y-axis. 
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(a)  Elevation  Angle  Jacobians  for  the  n- 
North,  East,  Down  axes  respectively 
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(b)  Scaled  View  of  Figure  4.35a 


Figure  4.35.  Features  at  Moon  Distance  along  E-axis  and  ECEF  Y-axis. 
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(a)  Elevation  Angle  Jacobians  for  the  n- frame 


(b)  Scaled  View  of  Figure  4.36a 


North,  East,  Down  axes  respectively 


Figure  4.36.  Features  at  Sun  Distance  along  E-axis  and  ECEF  Y-axis. 


n-frame  D-axis  Distance  Testing. 

For  this  set  of  scenarios,  lunar  features  were  placed  along  the  D-axis  at  great 
distances  to  observe  the  INS  error  correction  under  such  circumstances.  The  location 
of  the  lunar  features  along  the  ECEF  X-axis  for  this  set  of  scenarios  coincides  directly 
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with  the  n-frame  D-axis  as  the  initial  location  of  the  aircraft’s  trajectory  is  located 
at  0°lat  0°lon  and  stays  within  a  small  range  over  its  course  as  shown  in  Figure 
4.37.  This  feature  orientation  with  respect  to  the  aircraft  provides  for  a  special  case 
in  position  state  observabilities  provided  by  the  Jacobians  due  to  the  nature  of  the 
spherical  coordinate  system  angles.  The  azimuth  angle  Jacobians  are  solely  calculated 
for  the  north  and  east  position  states  with  observabilities  along  the  planar  surface  of 
the  E-N  plane  independent  of  the  down  position  state  observability.  The  E-N  plane  is 
perpendicular  to  the  D-axis,  essentially,  still  holding  Jacobians  with  observability  of 
the  north  and  east  position  states  independent  of  the  feature  distances  along  the  D- 
axis.  The  elevation  angle  Jacobians  are,  however,  tied  to  the  position  of  the  features 
along  the  D-axis  based  upon  the  nature  of  the  elevation  angle  holding  observabilties 
along  all  three  position  states. 
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Figure  4.37.  ECEF  Frame.  Features  along  n-frame  D-axis  and  ECEF  X-axis  (Not  to 
Scale). 

Figures  4.38a  and  4.39a  present  the  results  of  the  INS  drift  correction  based  upon 
lunar  features  directly  above  the  n-frame  along  the  D-axis  coinciding  with  the  ECEF 
X-axis.  Both  sets  of  error  variance  for  the  north  and  east  position  states  stay  within 
±0.35m  when  the  features  are  moved  from  the  lunar  distance  above  n-frame  to  the 
Sun’s  distance.  The  behavior  is  reiterated  in  the  ensemble  statistics  of  the  Jacobians. 
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The  north  and  east  position  state’s  azimuth  angle  Jacobians  with  features  at  lunar  and 
solar  distances  in  Figures  4.40  and  4.41  respectively  do  not  show  a  dramatic  change 
in  observability  with  the  change  in  distance  along  the  D-axis.  On  the  other  hand,  the 
position  state  observabilities  from  the  elevation  angle  Jacobians  in  Figures  4.42  and 
4.43  do  show  a  greater  decrease  by  an  order  of  magnitude  for  the  all  three  position 
states,  but  the  most  dramatic  decrease  is  seen  in  the  change  for  the  observability  of 
the  down  position  state  by  the  elevation  angle  Jacobian. 


(a)  Position  State  Errors  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.38.  Features  at  Moon  Distance  along  D-axis  and  ECEF  X-axis. 
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(a)  Position  State  Errors  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.39.  Features  at  Sun  Distance  along  D-axis  and  ECEF  X-axis. 
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(a)  Azimuth  Angle  Jacobians  for  the  n-frame 
North,  East,  Down  axes  respectively 
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(b)  Scaled  View  of  Figure  4.40a 


Figure  4.40.  Features  at  Moon  Distance  along  D-axis  and  ECEF  X-axis. 
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(a)  Azimuth  Angle  Jacobians  for  the  n- frame 
North,  East,  Down  axes  respectively 
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(b)  Scaled  View  of  Figure  4.41a 


Figure  4.41.  Features  at  Sun  Distance  along  D-axis  and  ECEF  X-axis. 
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(a)  Elevation  Angle  Jacobians  for  the  n-frame  (b)  Scaled  View  of  Figure  4.42a 

North,  East,  Down  axes  respectively 


Figure  4.42.  Features  at  Moon  Distance  along  D-axis  and  ECEF  X-axis. 
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(a)  Elevation  Angle  Jacobians  for  the  n-frame  (b)  Scaled  View  of  Figure  4.43a 

North,  East,  Down  axes  respectively 


Figure  4.43.  Features  at  Sun  Distance  along  D-axis  and  ECEF  X-axis. 


Lunar  Orbital  Geometry  Testing. 

The  following  set  of  scenarios  presents  the  results  of  performing  INS  drift  cor¬ 
rection  aided  by  lunar  feature  angle  measurements  under  different  lunar  geometries 
obtained  from  lunar  orbit  characteristics  detailed  in  Table  2.1.  The  previous  section 
showed  that  features  located  at  great  distances  along  an  n-frame  axis  will  cause  mini¬ 
mal  INS  drift  correction  for  the  associated  position  state  error.  With  this  in  mind,  the 
error  state  estimation  and  correction  performance  of  the  EKF  under  different  lunar 
orientations  can  be  predicted.  As  mentioned  in  Section  2.6,  the  two  main  categories 
of  lunar  celestial  geometry  that  will  be  looked  at  are: 

1.  Lunar  feature  location  at  inclination  with  respect  to  the  Earth’s  equator  plane. 

2.  Lunar  feature  location  along  the  lunar  orbit  path  with  respect  to  the  Earth’s 
center  and  Prime  Meridian. 
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Inclination  to  Equator. 


Different  lunar  inclinations  to  the  equator  are  shown  in  Figure  4.44  providing  the 
frame  for  the  following  scenarios: 

1.  Lunar  features  at  the  highest  inclination  to  equator. 

2.  Lunar  features  above  the  equator  and  at  the  zenith  of  the  trajectory. 

3.  Lunar  features  at  the  lowest  inclination  to  equator. 
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Figure  4.44.  Lunar  Inclinations  to  Equator. 


For  the  lunar  features  at  the  highest  inclination  to  the  equator,  Figures  4.45- 
4.47  show  the  orientation  of  the  simulated  lunar  features  with  respect  to  the  ECEF 
frame  and  the  n- frame.  Figure  4.48a  shows  a  greater  error  variance  in  the  north 
position  state  at  ±12. 13m  along  with  lower  variance  for  the  east  position  state  error 
at  ±3. 981m. 
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Figure  4.45.  ECEF  Frame.  Highest  Lunar  Inclination  to  Equator  (Not  to  Scale). 


n-frame 


X  Position  (km) 


Figure  4.46.  ECEF  Frame:  X-Z  Plane.  Highest  Lunar  Inclination  to  Equator  (Not  to 
Scale). 
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Y  Position  (km) 


Figure  4.47.  ECEF  Frame:  Y-Z  Plane.  Highest  Lunar  Inclination  to  Equator  (Not  to 
Scale). 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.48.  Features  at  Highest  Lunar  Inclination  to  Equator. 
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Figures  4.49  and  4.50  present  the  scenario  where  lunar  features  are  located  straight 
above  the  equator  and  the  trajectory  zenith.  The  position  state  errors  are  shown  in 
Figure  4.51a.  Both  the  east  and  north  position  error  states  are  within  ±0.33m.  The 
geometry  of  the  angles  straight  above  the  trajectory  provide  for  excellent  observability 
of  the  north  and  east  position  states  which  he  along  the  N-E  plane  with  the  features 
lying  directly  along  the  D-axis  above.  This  follows  the  special  case  of  the  spherical  co¬ 
ordinate  system  angle  calculations  mentioned  in  the  previous  n-frame  D-axis  distance 
testing  subsection. 


I 


Figure  4.49.  ECEF  Frame:  Inclination  to  Equator  at  0°  (Not  to  Scale). 


n-frame 


X  Position  (km) 


Figure  4.50.  ECEF  Frame:  X-Z  Plane.  Inclination  to  Equator  at  0°  (Not  to  Scale). 
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(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 


Figure  4.51.  Features  above  Equator  and  at  Trajectory  Zenith. 


The  next  scenario  involves  the  behavior  of  the  position  state  estimates  aided  by 
lunar  feature  measurements  located  at  the  lowest  inclination  to  the  equator  repre¬ 
sented  by  Figure  4.52.  The  results  can  be  seen  in  Figures  4.53a.  It  is  interesting 
to  note  that  the  north  and  east  position  state  error  variance  for  lunar  features  at 
the  highest  inclination  to  the  equator  in  Figure  4.48a  follows  a  similar  error  variance 
pattern  for  the  position  state  errors  for  the  lunar  features  at  the  lowest  inclination  to 
the  equator.  This  can  be  explained  by  both  scenario  sets  having  features  located  at 
+28.58°  or  —28. 58° along  the  same  plane  used  in  both  scenarios  allowing  for  similar 
Jacobians  to  be  calculated  for  the  position  states. 
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Figure  4.52.  ECEF  Frame:  Lowest  Inclination  to  Equator  (Not  to  Scale). 


(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 

Figure  4.53.  Features  at  Lowest  Lunar  Inclination  to  Equator. 


Lunar  Feature  Position  Along  the  Orbit  with  respect  to  the  Prime 
Meridian. 

The  following  scenarios  are  constructed  to  understand  the  behavior  of  the  INS  drift 
correction  aided  by  lunar  features  that  are  located  at  positions  along  the  lunar  orbit 
with  respect  to  the  Earth’s  center  and  Prime  Meridian  as  represented  in  Figure  4.54. 
The  first  scenario  is  represented  in  Figure  4.55  where  the  location  of  the  trajectory  is 
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at  (Flat  O°lon  and  the  lunar  features  lie  along  the  Moon’s  orbit  at  90°W  of  the  Prime 
Meridian  still  maintaining  the  default  high  inclination  to  the  equator  discussed  in 
Section  4.1.  Based  upon  the  behavior  of  the  EKF  and  the  measured  angles,  the 
estimated  east  position  state  will  have  a  higher  error  variance  than  the  north  position 
state  from  what  was  discovered  in  the  feature  distance  testing  with  regards  to  angle 
observability  of  position  states. 


Moon  above  90°W 
longitude 


Figure  4.54.  Lunar  feature  orientations  to  Earth  longitude. 
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The  results  from  the  scenario  can  be  viewed  in  Figures  4.56a.  As  expected,  the 
error  variance  for  the  east  position  state  is  much  larger  at  ±228. 2m  than  the  north 
position  state  at  ±123. 7m. 


(a)  Position  State  Error  for  the  n-frame  North, 
East,  Down  axes  respectively 


(b)  Estimated  vs.  True  Trajectory 


Figure  4.56.  Lunar  features  at  90°W  of  the  Prime  Meridian  and  Highest  Inclination  to 
Equator. 
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The  next  scenario  involved  the  lunar  features  being  placed  along  the  Moon’s  orbit 
at  45°W  of  the  Prime  Meridian  and  at  the  highest  inclination  to  the  equator  28.58°  as 
shown  in  Figure  4.57.  The  results  can  be  viewed  in  Figures  4.58a.  The  results  show 
a  decrease  in  error  variance  for  the  north  and  east  position  states  from  the  results 
in  the  previous  scenario  with  lunar  features  located  at  90°W.  Figure  4.58a  shows  the 
north  position  state  now  has  ±13. 97m  and  the  east  position  state  has  ±16. 71m  error 
variance. 


£  *io<  , 

n-frame 


Figure  4.57.  ECEF  Frame:  Lunar  features  at  45°W  and  Highest  Inclination  to  Equator 
(Not  to  Scale). 
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(a)  Position  State  Error  for  the  n-frame  North,  (b)  Estimated  vs.  True  Trajectory 

East,  Down  axes  respectively 

Figure  4.58.  Lunar  features  at  45°W  of  the  Prime  Meridian  and  Highest  Inclination  to 
Equator. 

The  scenario  for  the  features  above  the  Prime  Meridian  and  at  the  highest  in¬ 
clination  to  the  equator  holds  the  same  orientation  for  the  trajectory  and  features 
represented  in  Figures  4.45-4.47.  With  the  features  at  this  location,  the  error  variance 
for  the  north  position  state  decreased  to  ±12. 13m,  and  the  east  position  state  error 
variance  decreased  to  ±3. 981m  as  shown  in  Figure  4.48a. 

4.6  Summary 

Chapter  4  presented  the  results  of  the  research  methodology  defined  in  Chapter 
3.  The  performance  of  an  EKF  to  estimate  the  error  of  an  INS  with  the  aiding 
of  lunar  features  was  investigated  under  various  scenarios.  Section  4.1  laid  out  the 
default  trajectory  parameters  and  the  default  orientation  between  the  lunar  features, 
aircraft  trajectory,  and  Earth  in  the  ECEF  frame  and  n-frame.  Section  4.2  showed 
the  effect  of  aiding  by  lunar  feature  angle  measurements  on  different  INS  grades. 
Section  4.3  tested  INS  drift  correction  aided  by  lunar  feature  angle  measurements 
under  different  camera  specifications.  Section  4.4  looked  at  INS  drift  correction  aided 
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by  lunar  features  of  different  quantities.  The  effects  of  lunar  feature  distance  and 


geometry  on  lunar  vision  aided  navigation  with  INS  was  explored  in  Section  4.5. 
Section  4.6  provided  a  summary  for  the  chapter. 
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V.  Conclusions 


5.1  Introduction 

This  research  effort  presented  the  theoretical  limits  of  a  non-GPS  precision  navi¬ 
gation  method  consisting  of  an  INS  aided  by  lunar  feature  angle  measurements  where 
only  measurement  noise  effects  were  considered.  The  following  sections  summarize 
the  research  methods  and  results  of  the  scenarios  from  Chapters  3  and  4  respectively 
and  conclude  with  suggested  future  work. 

5.2  Conclusions 

The  purpose  of  this  research  effort  was  to  establish  the  feasibility  of  a  precision 
navigation  method  by  which  INS  drift  error  could  be  corrected  through  the  use  of 
measured  lunar  feature  angles.  To  accomplish  this  task,  an  EKF  algorithm  was  imple¬ 
mented  to  process  an  INS  trajectory  and  estimate  the  INS  error  states.  Corrections 
to  the  INS  error  states  were  made  during  the  EKF  update  stage  with  lunar  feature 
angle  measurements.  Through  the  use  of  Jacobians,  the  lunar  feature  angle  mea¬ 
surements  provided  position  state  observability  based  on  the  change  in  the  measured 
lunar  feature  angles  with  respect  to  the  change  in  distance  along  the  position  states. 
A  barometer  was  simulated  to  aid  the  down  position  error  state  with  an  update  rate 
of  60  seconds.  Four  total  categories  of  scenarios  were  accomplished  to  observe  the 
behavior  of  the  INS  drift  corrections  under  different  lunar  feature  conditions. 

The  first  scenario  set  showed  marked  improvement  in  the  navigation  solution 
error  at  la  for  an  INS  aided  by  lunar  feature  angle  measurements  compared  to  a 
stand  alone  INS.  Also,  the  scenario  set  showed  that  the  increase  in  INS  grade  quality 
led  to  improved  position  state  estimation.  Table  5.1  presents  the  position  state  error 
results  for  the  first  scenario  set.  The  use  of  Commercial  Grade  INS  in  the  EKF  with 


103 


100  lunar  feature  angle  measurements  brought  about  position  state  error  variance  of 
±58. 96m  for  the  north  position  state  and  ±25. 85m  for  the  east  position  state.  When 
using  the  Navigation  Grade  INS,  the  position  state  error  decreased  significantly  to 
±2. 946m  for  the  north  position  state  error  variance  and  ±0.9805  for  the  east  position 
state  error  variance. 

Table  5.1.  Feasibility  Table  for  Lunar- Aided  Navigation  with  respect  to  INS  grade. 
Navigation  Solution  Error  at  lcr  listed  in  the  right  three  columns  were  collected  at 
lOmin.  Parameters  are  with  respect  to  a  trajectory  traveled  in  a  circle  with  a  radius 
20km  radius  with  starting  location  at  0°lat  0°lon. 


INS 

Camera  Noise  (mrads) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

Commercial 

58.96 

25.85 

77.62 

Tactical 

0.011 

253 

100 

384,400 

Highest  Inclination  to  Equator 

13.5 

4.914 

16.43 

Navigation 

2.946 

0.9805 

2.328 

The  second  set  of  scenarios  looked  at  measurement  noise  effects  from  camera 
specifications  upon  INS  drift  corrections  by  lunar  feature  angle  measurements.  A 
sweep  of  the  measurement  noise  vopt^  based  on  different  camera  image  sensor  pixel 
sizes  A£  was  conducted.  As  shown  in  Table  5.2,  at  the  largest  pixel  size  1000/xm 
which  equated  to  a  measurement  noise  of  2mrads,  the  error  variance  for  the  north 
position  state  was  ±396m  and  ±216. 5m  for  the  east  position  state.  Decreasing  the 
pixel  size  to  1/xm  and  the  associated  noise  to  0.002mrads  caused  the  INS  error  solution 
to  improve  to  an  error  variance  of  ±8. 33m  for  the  north  position  state  and  ±  1.54m 
for  the  east  position  state. 

Table  5.2.  Feasibility  Table  for  Lunar- Aided  Navigation  with  INS  with  respect  to 
camera  image  sensor  pixel  size.  Navigation  Solution  Error  at  la  listed  in  the  right 
three  columns  were  collected  at  lOmin.  Parameters  are  with  respect  to  a  trajectory 
traveled  in  a  circle  with  a  radius  20km  radius  with  starting  location  at  0°lat  0°lon. 


INS 

Pixel  Size  (/im) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

Tactical 

1000 

253 

100 

384,400 

Highest  Inclination  to  Equator 

396 

216.5 

13.18 

100 

54.04 

31 

12.53 

10 

16.97 

6.473 

12.68 

1 

8.33 

1.54 

14.43 
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Parameters  from  a  commercial  camera  were  obtained  to  see  the  effects  of  different 


lens  focal  lengths  on  the  measurement  noise  and  position  state  error  corrected  by 
lunar  feature  angle  measurements.  Overall,  the  results  in  Table  5.3  pointed  to  the 
increase  in  lens  focal  length  decreasing  the  position  state  error  and  increasing  the 
performance  of  the  EKF  in  estimating  the  true  trajectory. 

Table  5.3.  Feasibility  Table  for  Lunar- Aided  Navigation  with  INS  with  respect  to  Lens 
Focal  Length.  Navigation  Solution  Error  at  la  listed  in  the  right  three  columns  were 
collected  at  lOmin.  Parameters  are  with  respect  to  a  trajectory  traveled  in  a  circle 
with  a  radius  20km  radius  with  starting  location  at  0°lat  0°lon. 


INS 

Camera  Noise  (mrads) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

0.055 

50 

28.67 

11.25 

10.55 

Tactical 

0.011 

253 

100 

384,400 

Highest  Inclination  to  Equator 

7.62 

4.805 

11.36 

0.003 

1000 

6.42 

2.347 

9.88 

The  third  scenario  set  showed  that  an  increase  in  the  lunar  feature  quantity  de¬ 
creased  position  state  error  resulting  in  a  stronger  trajectory  estimate.  The  results  of 
the  third  scenario  set  are  presented  in  Table  5.4.  At  one  measured  lunar  feature,  the 
error  variance  was  at  ±273. 9m  for  the  north  position  state  and  ±114m  for  the  east 
position  state.  At  1000  features,  the  INS  solution  was  corrected  to  an  error  variance 
of  ±7. 435m  for  the  north  position  state  and  ±1.326m  for  the  east  position  state. 

Table  5.4.  Feasibility  Table  for  Lunar- Aided  Navigation  with  INS  with  respect  to 
Feature  Quantity.  Navigation  Solution  Error  at  la  listed  in  the  right  three  columns 
were  collected  at  lOmin.  Parameters  are  with  respect  to  a  trajectory  traveled  in  a 
circle  with  a  radius  20km  radius  with  starting  location  at  0°lat  0°lon. 


INS 

Camera  Noise  (mrads) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

Tactical 

0.011 

253 

1 

384,400 

Highest  Inclination  to  Equator 

273.9 

114 

14.56 

10 

33.09 

26.61 

15.69 

100 

13.5 

4.914 

16.43 

1000 

7.435 

1.326 

13.31 

The  next  set  of  scenarios  encompassed  the  effects  of  feature  distance  and  different 
lunar  orbit  geometries  on  INS  drift  correction.  In  Section  4.5,  the  feature  distance 
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scenario  set  showed  that  features  at  distances  approaching  infinity  will  cause  a  lack  of 
observabilty  of  the  position  states  resulting  in  an  increased  position  state  error.  Also, 
the  feature  distance  set  of  scenarios  in  Section  4.5  confirmed  the  assertion  that  a  lack 
of  observability  will  occur  for  a  distinct  position  state  when  features  lie  at  far  distances 
along  the  respective  position  axis.  However,  for  the  increased  feature  distance  along 
the  n-frame  D-axis  at  the  trajectory’s  zenith,  minimal  effect  occurred  to  the  position 
state  errors.  Also,  at  this  orientation,  the  north  and  east  position  state  error  variances 
were  both  shown  to  be  extremely  small  ±0.35m.  This  was  due  in  part  to  the  azimuth 
Jacobians  providing  observabilities  of  the  north  and  east  position  states  independent 
of  the  down  position  state.  The  error  results  for  the  distance  scenarios  are  not  listed 
in  a  table  as  the  feature  locations  are  not  set  at  realistic  orbit  parameters.  The  error 
results  from  the  distance  scenarios  can  be  found  in  Section  4.5. 

For  the  scenario  set  dealing  with  INS  drift  correction  based  on  lunar  features 
at  different  orbital  geometries,  features  were  placed  at  different  inclinations  to  the 
equator  and  at  different  locations  along  the  lunar  orbit  with  the  respect  to  the  Prime 
Meridian.  The  lunar  geometries  are  based  on  orbital  parameters  listed  in  Table  2.1. 
During  the  changes  in  feature  locations  at  different  lunar  geometries,  the  location  of 
the  true  trajectory  maintained  a  circular  path  with  a  20  km  radius  holding  a  starting 
position  at  0°lat  0°lon.  Table  5.5  shows  the  position  state  error  results  for  the  lunar 
surface  features  placed  at  different  lunar  inclinations  with  respect  to  the  Earth’s 
equator.  The  position  state  error  for  the  north  and  east  position  states  were  shown 
to  be  minimally  changed  when  the  inclination  of  the  feature  locations  at  28.58°  was 
mirrored  and  the  features  were  placed  on  the  Moon’s  orbit  at  the  lowest  inclination 
to  the  equator.  This  is  due  to  the  similar  change  in  angles  occurring  at  both  feature 
locations  causing  similar  Jacobians  which  provide  the  observability  of  the  states  to 
aid  in  error  state  estimation  correction. 
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Table  5.5.  Feasibility  Table  for  Lunar- Aided  Navigation  with  INS  with  respect  to  lunar 
surface  features  placed  at  different  Lunar  Inclinations  to  Earth’s  equator.  Navigation 
Solution  Error  at  la  listed  in  the  right  three  columns  were  collected  at  lOmin.  Param¬ 
eters  are  with  respect  to  a  trajectory  traveled  in  a  circle  with  a  radius  20km  radius 
with  starting  location  at  0°lat  0°lon. 


INS 

Camera  Noise  (mrads) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

Highest  Inclination  to  Equator 

12.13 

3.981 

11.88 

Tactical 

0.011 

253 

100 

384,400 

At  Zenith  &  Over  Equator 

0.322 

0.3275 

14.64 

Lowest  Inclination  to  Equator 

8.104 

4.098 

11.4 

With  features  placed  at  different  locations  along  the  lunar  orbit  with  respect  to 

the  Prime  Meridian  and  maintaining  the  highest  inclination  to  the  equator  at  28.58°, 

the  INS  error  correction  was  shown  in  Table  5.6  to  improve  as  the  features  approached 

the  zenith  of  the  trajectory.  With  the  lunar  features  located  at  90°W  of  the  Prime 

Meridian,  the  north  position  state  error  variance  was  ±123. 7m  and  the  east  position 

state  error  variance  was  ±228. 2m.  With  the  change  of  the  feature  locations  to  45°W 

of  the  Prime  Meridian,  the  north  and  east  position  state  error  variances  decreased  to 

±13. 97m  and  ±16. 71m  respectively.  The  decrease  in  error  continued  as  the  feature 

locations  were  moved  to  the  lunar  orbit  position  lying  above  the  Prime  Meridian 

with  an  inclination  of  28.58°  to  the  equator.  The  north  position  error  state  variance 

decreased  to  ±12. 13m  and  the  east  position  state  variance  decreased  to  ±3. 981m. 

Table  5.6.  Feasibility  Table  for  Lunar- Aided  Navigation  with  INS  with  respect  to 
lunar  surface  features  placed  along  different  lunar  orbit  positions  and  at  the  highest 
inclination  to  the  Equator’s  equator.  Navigation  Solution  Error  at  la  listed  in  the  right 
three  columns  were  collected  at  lOmin.  Parameters  are  with  respect  to  a  trajectory 
traveled  in  a  circle  with  a  radius  20km  radius  with  starting  location  at  0°lat  0°lon. 


INS 

Camera  Noise  (mrads) 

Focal  Length  (mm) 

#  of  features 

Distance  (km) 

Lunar  Geometry 

Error  Variance  (m) 

N 

E 

D 

90°W  of  Prime  Meridian 

123.7 

228.2 

6.02 

Tactical 

0.011 

253 

100 

384,400 

45°W  of  Prime  Meridian 

13.97 

16.71 

11.37 

Over  Prime  Meridian 

12.13 

3.981 

11.88 

In  conclusion,  this  research  effort  presented  the  theoretical  limits  of  a  viable  non- 
GPS  precision  navigation  method  able  to  operate  when  GPS  is  unavailable.  The 
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method  is  based  on  lunar  feature  angle  measurements  aiding  an  INS.  Although  mea¬ 
surement  noise  was  only  considered  to  effect  the  model,  the  results  in  Chapter  4 
showed  that  INS  drift  correction  provided  by  lunar  feature  angle  measurements  can 
theoretically  greatly  improve  a  navigation  solution  compared  to  a  stand-alone  INS. 

5.3  Future  Work 

With  the  conclusion  of  this  research  effort,  below  are  areas  that  are  worth  exploring 
in  continuing  to  pursue  lunar  vision-aided  navigation: 

1.  Conduct  INS  drift  correction  with  lunar  feature  angle  measurements  under  re¬ 
alistic  noise  parameters  considering  error  in  feature  extraction  and  matching 
techniques,  and  the  error  for  a  camera’s  alignment  to  the  local  navigation  frame. 

2.  Observe  the  simulated  INS  drift  correction  aided  by  feature  angles  extracted 
from  lunar  feature  images.  Capture  real  world  lunar  images  with  a  camera  to 
understand  the  relationship  between  camera  resolution  and  the  observed  lunar 
images. 

3.  Implement  a  feature  detection  and  extraction  algorithm  such  as  SIFT  or  SURF 
on  lunar  images  to  find  out  which  lunar  surface  characteristics  (seas,  craters, 
etc.)  present  the  most  robust  source  of  features  to  use  for  angle  measurements. 
Investigate  feature  detection  and  extraction  performance  under  different  visi¬ 
bility  scenarios  such  as  moon  phases,  occlusion,  librations,  etc. 
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